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LEXINGTON 


MASSACHUSETTS 


ABSTRACT 


This  report  describes  25  investigations  in  the  field  of  seismic 
discrimination.  These  are  grouped  as  follows:  estimation  of 

body  wave  magnitude  by  single  stations  and  networks  (9  con- 
tributions), surface  wave  studies  (2  contributions),  examination 
of  methods  for  the  separation  of  the  signals  from  a sequence  of 
closely  spaced  explosions  (2  contributions),  studies  concerned 
with  lateral  heterogeneity  within  the  earth  (5  contributions), 
miscellaneous  studies  (3  contributions),  and  recent  developments 
in  our  data  and  computer  systems  (4  contributions). 


CONTENTS 


Abstract  iii 

Summary  vii 

I.  BODY  WAVE  MAGNITUDE  mb  1 

A.  Global  Variations  in  Station  Magnitude  Bias  1 

B.  Source  Region  Variations  in  Bias  4 

C.  Temporal  Variations  in  the  Shape  of  the  Magnitude- 

Frequency  Curve  5 

D.  Comparison  of  Recent  Estimates  of  Magnitude  Bias  7 

E.  Seismic  Magnitude  Simulation  Studies  9 

1.  True  Magnitude  Distributions  9 

2.  Bias  and  Scatter  Parameters  10 

3.  Station  Detection  Characteristics  10 

F.  Magnitude  Differences  Between  Station  Pairs  14 

G.  Estimation  of  Seismic  Source  at  1 Hz  from  log(A/T)  15 

H.  Evidence  for  a; -Cube  Scaling  from  Short-Period  Amplitudes 

and  Periods  of  the  Rat  Island  Sequence  (1965)  17 

I . Contribution  of  WWSSN  Stations  to  the  PDE  Bulletin  19 

II.  SURFACE  WAVE  STUDIES  39 

A.  The  Effects  of  Lateral  Heterogeneity  in  the  Oceanic 

Lithosphere  Upon  Surface  Wave  Propagation  39 

B.  Synthetic  Rayleigh  Wave  Seismograms  for  the  Fallon 

Earthquake  40 

III.  EVASION  49 

A.  Detecting  Multiple  Events  in  Narrowband  Seismograms  49 

B.  Comparison  of  Different  Prediction  Error  Filters 

on  a Synthetic  Composite  Event  50 

IV.  EARTH  HETEROGENEITY  55 

A.  Analysis  of  Scattering  from  Novaya  Zemlya  Explosions  55 

B.  Mislocation  Patterns  for  Seismic  Networks  56 

C.  Large-Scale  Heterogeneities  in  the  Lower  Mantle; 

Correlation  with  the  Gravity  Field  57 

D.  Internal  Discontinuities  in  the  Earth's  Structure 

and  Perturbation  Theory  61 

E.  Transmission  Holography  with  Impulsive  Sources  6 5 


v 


V.  MISCELLANEOUS  STUDIES  79 

A.  Model  Order  Criteria  Applied  to  Autoregressive  and  Non- 

Autoregressive  Processes  79 

B.  Application  of  the  Energy- Moment  Tensor  Relation 

to  Finite  Duration  Sources  81 

C.  Contribution  of  Large  Earthquakes  to  Excitation  of  Polar 

Motion  for  Years  1901  to  1970  83 

VI.  DATA  SYSTEMS  9 5 

A.  Seismic  Data  and  Information  System  9 5 

B.  Multi -Azimuth  Large  Array  Long- Period  Data  Base  98 

C.  SRO  Data  Tape  Processing  Facility  98 

D.  PDP-11  Software  99 

Glossary  103 


vi 


SUMMARY 


This  is  the  twenty -fifth  Semiannual  Technical  Summary  report  describing  the  activities  of 
Lincoln  Laboratory,  M.I.T.,  in  the  field  of  seismic  discrimination.  These  activities  involve 
research  into  the  fundamental  seismological  problems  associated  with  the  detection,  location, 
and  identification  of  earthquakes  and  nuclear  explosions,  and  with  the  estimation  of  explosion 
yields.  In  order  to  investigate  these  problems,  we  are  continuously  improving  our  ability  to 
manipulate  and  process  seismic  data  from  a global  network  of  high-quality  digital  stations  and 
arrays. 

We  are  engaged  in  a series  of  related  projects,  all  aimed  at  improving  our  understanding 
of  body  wave  magnitude  m^,  a quantity  of  considerable  importance  for  teleseismic  yield  esti- 
mation. In  particular,  we  have  focused  our  attention  on  the  various  kinds  of  bias  that  can  affect 
single-station  and  network  measurements  of  m^  for  an  event.  One  study  has  reduced  the  large 
volume  of  magnitude  data  contained  in  the  ISC  Bulletin,  in  order  to  obtain  values  of  station  bias 
for  a set  of  72  stations.  The  resulting  biases  show  good  correlation  with  tectonic  region,  partic- 
ularly in  the  U.S.  At  the  same  time,  the  results  of  this  study  indicate  disturbing  nonuniformity 
in  operator  performance  at  a number  of  stations. 

Earlier  work  on  attempting  to  remove  bias  due  to  detection  characteristics  is  continuing. 
We  have  begun  a series  of  numerical  simulations  which  will  allow  us  to  examine  the  contributions 
of  the  various  different  parameters,  including  seismicity  models  and  instrument  saturation.  At 
the  same  time,  we  are  developing  various  methods  for  the  estimation  of  these  parameters.  The 
use  of  magnitudes  measured  with  respect  to  a reference  station  appears  particularly  promising. 
Related  studies  are  concerned  with  removing  the  bias  that  can  arise  from  a combination  of  the 
shape  of  the  source  spectrum,  attenuation  in  the  earth,  and  instrument  response.  The  tradi- 
tional method  of  dividing  the  observed  amplitude  by  the  dominant  period  of  the  signal  is  only  a 
naive  approach  to  a complex  problem.  Indications  are  given  that  it  will  be  possible  to  make  a 
much  more  sophisticated  approach  to  the  removal  of  these  effects. 

Improvements  to  previous  Lincoln  Laboratory  work  on  the  use  of  ray  tracing  techniques  to 
study  the  propagation  of  surface  waves  from  sources  in  Eurasia  have  been  made.  Lateral  heter- 
ogeneity within  the  oceanic  lithosphere  has  been  included  by  utilizing  the  known  variation  of  phase 
velocity  with  age.  Other  surface  wave  studies  have  included  an  application  of  the  seismic  mo- 
ment tensor  formulism  to  the  analysis  of  earthquake  source  mechanism.  As  an  example,  the 
case  of  the  surface  waves  observed  at  LRSM  sites  from  the  Fallon  earthquake  is  described. 

We  are  concentrating  our  efforts  in  the  area  of  evasion  techniques  on  the  problem  of  the 
separation  of  the  signals  from  a closely  spaced  sequence  of  shots.  The  utilization  of  maximum 
entropy  cepstral  analysis  appears,  on  the  basis  of  a simulated  example,  to  be  promising.  In 
another  case,  adaptive  deconvolution  methods  are  quite  successful. 

A variety  of  studies  has  continued  our  efforts  to  characterize  lateral  heterogeneity  within 
the  earth.  Complexities  in  the  P codas  from  some  Novaya  Zemlya  explosions  are  shown  to 
originate  in  a broad  region  which  may  not  be  associated  with  source  or  receiver  region.  Earlier 
results  from  a global  analysis  of  travel-time  variations  have  been  extended.  It  is  shown  that  it 
is  possible  to  make  a convincing  correlation  between  the  inferred  velocity  anomalies  and  the 
gravitational  field  of  the  earth.  And  another  study  examines  the  determination  of  levels  of  ve- 
locity discontinuity  within  the  earth  using  free  oscillation  data. 


Vll 


Among  other  studies  reported  here  are  included  some  comparisons  between  various  criteria 
that  have  been  proposed  for  the  determination  of  the  optimum  prediction  error  filter  length  for 
maximum  entropy  spectral  analysis.  Also,  an  extension  of  previous  work  on  the  energy-moment 
tensor  relation  is  described.  In  another  investigation,  the  observed  time  series  for  the  excita- 
tion of  the  Chandler  motion  of  the  earth  is  compared  with  the  calculated  contribution  due  to  large 
earthquakes.  A good  correlation  is  observed. 

We  have  continued  to  develop  our  data  and  computer  systems  in  order  to  be  able  to  handle 
the  SRO  data  currently  on  hand,  and  the  much  larger  quantity  of  digital  data  soon  to  be  available. 
We  have  designed  an  expansion  of  our  CONSOLE  program  that  permits  interactive  display  and 
processing  of  SRO  data  on  our  PDP-7  computers.  Our  PDP-11  based  system  is  already  operat- 
ing as  an  ARPANET  connection,  and  will  become  a general-purpose  time -sharing  facility  in  the 
very  near  future. 


M.  A.  Chinnery 


SEISMIC  DISCRIMINATION 


I.  BODY  WAVE  MAGNITUDE  mb 

A.  GLOBAL  VARIATIONS  IN  STATION  MAGNITUDE  BIAS 
1 

In  the  last  SATS,  the  results  of  a study  of  body  wave  magnitudes  reported  in  the  USCGS 
earthquake  data  reports  (EDR)  lists  for  the  first  6 months  of  1971  were  given.  This  has  been 
extended  to  the  entire  set  of  station  m^  observations  given  in  the  ISC  Bulletins  for  1964-1973, 
comprising  over  400,000  m^  reports  for  nearly  40,000  events.  Mean  station  magnitude  biases 
have  been  computed  for  the  best  (in  terms  of  events  reported)  72  stations  with  respect  to  the 
mean  magnitude  of  observations  reported  by  this  set  of  stations,  in  the  manner  previously  re- 
ported. The  only  difference  is  that  we  have  required  that  an  event  be  reported  by  more  than 
15  of  these  72  stations  before  a bias  is  calculated.  In  all  eases,  the  distribution  of  biases  is 
remarkably  close  to  normal  and  is,  thus,  well  characterized  by  its  mean  and  standard  deviation. 
In  Table  1-1,  the  size,  mean,  and  standard  deviation  of  the  bias  distribution  for  each  of  these 
72  stations  are  given.  The  mean  biases  range  from  40.37  (EDM,  Canada)  to  —0.32  m^  units  (TFO. 
Arizona)  and  the  standard  deviations  of  the  distributions  from  0.24  to  0.50.  The  means  them- 
selves are  very  well  determined  because  of  the  large  sample  sizes,  and  the  standard  deviation  of 
the  mean  in  no  case  exceeds  0.02  m^  units.  There  is,  thus,  no  doubt  that  the  large  \ ariations 
in  mean  bias  are  statistically  significant. 

Biases  also  have  been  calculated  for  each  individual  year  of  data  as  well  as  the  entire  time 
period  1964-1973.  In  general,  the  biases  show  little  temporal  variation,  but  there  are  some 
stations  for  which  the  variation  with  time  is  alarming.  EUR  (Eureka,  Nevada)  shows  an  almost 
linear  decrease  in  mean  bias  from  +0.3  to  —0.2  m^  units  over  1964-1973;  nearby  stations  in 
Utah  show  no  such  variation.  Substantial  changes  in  mean  bias  with  time  also  were  found  for 
arrays  in  Montana  (LASA)  and  Sweden  (Hagsfors);  the  years  of  largest  change  also  are  those  ot 
largest  variance  in  the  distribution  of  biases  and  it  appears  that  the  magnitude  reporting  of  these 
arrays  was  grossly  degraded  in  certain  years.  In  most  cases,  the  temporal  variations  in  bias 
are  incapable  of  rational  explanation. 

Biases  also  have  been  calculated  for  a further  30  stations  which  reported  fairly  frequently. 
The  sample  sizes  here  are  smaller  and  they  must  necessarily  be  considered  less  accurate;  in 
particular,  we  cannot  determine  whether  they  vary  with  time  or  with  source  region.  The  size, 
mean,  and  standard  deviation  of  bias  for  these  30  stations  are  given  in  Table  1-2. 

Figure  1-1  shows  the  mean  biases  measured  for  the  total  of  102  stations  as  a function  of 
geographical  position.  Negative  biases  imply  low  Q,  and  positive  biases  high  Q,  beneath  the 

station.  It  can  be  seen  that  the  biases  are  in  accordance  with  other  measurements  of  Q from, 

2 3 

e.g.,  S propagation,  long-period  teleseismic  P-  and  S-wave  amplitudes,  and  surfaee  wave 

n 4 

attenuation.  Q is  lowest  in  the  region  of  mid-ocean  ridges  and  ’’rift"  structures  such  as  the 
western  U.S.  and  East  Africa,  and  highest  in  stable  regions  such  as  shields,  aseismic  platforms, 
and  deep  ocean  basins. 

Figure  1-2  shows  the  biases  for  stations  in  the  continental  U.S.  on  a larger  scale.  The  dif- 
ference in  attenuation  between  the  western  and  eastern  U.S.  is  immediately  apparent.  The 
dashed  lines  separate  regions  of  negative  and  positive  bias.  Notable  exceptions  to  this  rule  are 
the  Vela  stations  WMO  and  CPO;  this  may  be  partly  due  to  the  response  characteristics  of  the 
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TABLE  1-1 

MEAN  BIASES  FOR  72  STATIONS 

Number  of  Eventst 

Standord 

Deviation^ 

Number  of  Eventst 

Stondord 

Deviotiont 

Station 

N3 

N,5 

Mean  Bias 

Station 

N3 

N15 

Mean  Bios 

ALE 

1390 

1160 

-0.04 

0.29 

KRR 

3159 

1047 

-0. 24 

0.30 

ALQ 

5206 

1917 

-0.20 

0. 33 

KTG 

2221 

1798 

+ 0.  02 

0.30 

ASP 

1813 

456 

-0.  05 

0.  35 

LAO 

10202 

1705 

-0.  10 

0.47 

BHA 

2829 

980 

-0.28 

0.32 

LJU 

2501 

1739 

+ 0.29 

0.30 

BKS 

1476 

877 

+ 0.  18 

0.38 

LON 

2093 

1333 

-0.30 

0.37 

BMO 

26192 

3777 

-0.  29 

0.35 

LOR 

1879 

1147 

+ 0.06 

0.42 

BNG 

3944 

1029 

-0.07 

0.50 

LPS 

1448 

679 

+ 0.04 

0. 34 

BNS 

1975 

1695 

+ 0.20 

0.29 

MBC 

2327 

1256 

+ 0.  14 

0.  34 

BOZ 

2748 

978 

-0.06 

0.31 

MOX 

4645 

2762 

+ 0.  02 

0.27 

BUL 

4793 

1368 

-0. 07 

0.29 

MUN 

1738 

397 

+ 0.15 

0.37 

CAN 

2462 

507 

-0.02 

0.40 

NAO 

3399 

755 

-0.09 

0. 29 

CAR 

1289 

674 

+ 0.13 

0.38 

NDI 

1653 

906 

+ 0.33 

0.37 

CIR 

2623 

938 

-0.27 

0.30 

NEW 

2116 

1221 

+ 0.05 

0.30 

CLL 

2726 

1715 

+ 0.20 

0.32 

N1E 

1566 

1096 

-0.02 

0.33 

CLK 

2451 

902 

-0. 27 

0. 28 

NOR 

3693 

1986 

-0. 14 

0.33 

COL 

11418 

2579 

+ 0.01 

0.33 

NP- 

3934 

1035 

-0.00 

0.38 

COP 

1457 

1290 

+ 0. 36 

0.26 

NUR 

7256 

3149 

+ 0. 19 

0.30 

CPO 

11569 

2510 

-0.  07 

0. 35 

PMG 

4435 

1146 

+ 0. 10 

0.38 

DUG 

6117 

1753 

-0. 15 

0. 35 

PMR 

6016 

2075 

-0.08 

0.37 

EDM 

1709 

1181 

+ 0.37 

0.  28 

PNT 

1178 

663 

+ 0.13 

0.30 

EKA 

3094 

1395 

+ 0.00 

0.33 

POO 

2841 

1288 

+ 0. 17 

0.36 

EUR 

14997 

2913 

-0. 24 

0. 40 

PRE 

1947 

802 

-0. 07 

0. 39 

FUR 

2498 

1816 

+ 0.10 

0.38 

PRU 

2779 

2210 

+ 0.04 

0.  24 

GDH 

1439 

1257 

+ 0.00 

0.34 

RES 

1435 

1064 

+ 0.  13 

0.37 

GIL 

6138 

1842 

-0.04 

0. 35 

RIV 

1730 

507 

+ 0.31 

0.33 

GOL 

4058 

1770 

-0.  28 

0.39 

SHL 

2753 

769 

+ 0.  11 

0. 33 

GRF 

1450 

1104 

+ 0.  24 

0.28 

SJG 

1547 

718 

+ 0. 24 

0.38 

HFS 

2917 

895 

+ 0.05 

0. 45 

STU 

1752 

1434 

+ 0.29 

0.31 

HYB 

2106 

1166 

+ 0.  19 

0.37 

TFO 

18546 

2395 

-0.32 

0.  35 

KEV 

3678 

2218 

+ 0.  02 

0.27 

TRN 

1337 

704 

+ 0. 07 

0. 35 

KHC 

4092 

2720 

+ 0.  10 

0. 26 

TSK 

2762 

1047 

-0.  07 

0. 45 

KJF 

2137 

1036 

+ 0.09 

0.28 

TUC 

3033 

1263 

-0.14 

0.  25 

KJN 

6940 

2400 

+ 0.  14 

0.30 

TUL 

2820 

1152 

+ 0.21 

0.32 

KOD 

2294 

1210 

+ 0.06 

0.31 

UBO 

23157 

2828 

-0.11 

0.38 

KON 

1442 

1093 

+ 0.  07 

0.30 

WIN 

1259 

643 

-0.09 

0. 43 

KRA 

1743 

1251 

+ 0. 22 

0.29 

WMO 

13101 

1658 

-0.17 

0.31 

t Ng:  number  of  events  reported  by  3 or  more  stations. 
N^:  number  of  events  reported  by  15  or  more  stotions. 

tStondard  deviotion  of  distribution  of  magnitude  bioses. 

2 


TABLE  1-2 


MEAN  BIASES  FOR  30  SUPPLEMENTARY  STATIONS 


Station 

Number 
af  Eventst 

Mean  Bias 

Standard 

Deviatian-t 

ABU 

702 

+ 0.29 

0.53 

ATL 

213 

+ 0. 35 

0.51 

BRG 

749 

-0.  11 

0.26 

BUH 

414 

-0.04 

0.38 

DAG 

394 

-0.  02 

0.27 

ESK 

650 

+ 0.  19 

0.37 

FFC 

568 

+ 0.  08 

0.29 

FLO 

495 

+ 0. 26 

0.31 

GBA 

350 

+ 0.04 

0.37 

ILG 

314 

+ 0.04 

0.  42 

INK 

904 

+ 0.  15 

0.29 

KBL 

578 

+ 0.09 

0.31 

KLG 

639 

-0.02 

0.45 

LHN 

476 

-0. 15 

0.  36 

LPB 

455 

+ 0.  07 

0.31 

MAW 

190 

+ 0.  11 

0.  35 

MTD 

432 

-0.22 

0.30 

OIS 

392 

-0.  08 

0.41 

PAE 

248 

-0.03 

0.29 

PMO 

277 

+ 0.08 

0.35 

PNS 

370 

-0. 08 

0.45 

PPN 

233 

-0.  19 

0.32 

PPT 

237 

+ 0.06 

0.  29 

RAB 

766 

+ 0. 12 

0.43 

RCD 

264 

+ 0.  24 

0.31 

ROL 

392 

1 +0.19 

0.31 

RUV 

227 

+ 0.09 

0.37 

TPT 

299 

+ 0.06 

0.34 

TVO 

235 

+ 0.  08 

0.  28 

VAH 

283 

+ 0.  02 

0.  32 

t Number  of  events  reported  by  15  ar  mare  stations. 
t Standard  deviation  of  distribution  af  magnitude  biases. 
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instruments,  which  are  capable  of  recording  much  shorter  period  signals  than  conventional 
(e.g.,  WWSSN)  equipment.  Division  of  reported  amplitude  A by  dominant  period  T does  not 
entirely  compensate  for  the  much  higher  attenuation  at  these  shorter  periods. 

An  obvious  application  of  these  biases  is  their  use  as  station  magnitude  corrections  to  re- 
duce scatter  in  magnitude  observations  for  a single  event.  We  define  the  "scatter"  of  station 
m^  observations  relative  to  the  average  m^  by 
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for  events  in  a particular  magnitude  range  + Am^,  where  is  the  total  number  of 

events  in  the  magnitude  range  with  a total  of  associated  station  observations,  there  being  N. 
(>3)  observations  m.j  for  each  event  of  average  magnitude  m-.  Figure  1-3  shows  the  scatter 
as  defined  above,  with  Am^  = 0.1  units,  for  the  38,316  events  of  R > 3 during  1964-1973, 
prior  to  and  after  application  of  the  mean  station  biases  as  magnitude  corrections.  The  biases 
have  been  obtained  from  the  4668  events  with  IvL  >15,  nearly  all  5.0  < < 6.5.  Application 

of  the  corrections  has  decreased  the  scatter  by  at  least  15  percent  not  only  in  this  magnitude 
range  but  also  down  to  = 4.3.  Below  = 4.0,  the  scatter  has  been  increased;  this  is  pre- 
sumably because,  through  the  efforts  of  detection  characteristics,  our  biases  are  not  valid  at 
these  low  magnitudes.  Clearly,  the  scatter  could  be  further  reduced  by  applying  different  bias 
corrections  for  each  seismic  region  (and  year.1  ).  Unfortunately,  we  could  find  little  rational 
explanation  for  the  yearly  variation  in  station  bias,  and  no  consistency  in  regional  biases  across 
small  receiver  areas;  these  problems  should  be  resolved  before  more  sophisticated  station  cor- 
rections  are  applied. 


B.  SOURCE  REGION  VARIATIONS  IN  BIAS 


We  may  consider  the  mean  station  bias  b as  defined  in  Sec.  A,  above,  to  consist  of  the  fol- 
lowing factors 


b = bS  + bcs 


+ bUS  + bM  + b 


+ b^ 


c s cr 

where  b and  b are  introduced  by  crustal  structure  at  the  source  and  receiver,  respectively, 
us  u r 

b and  b by  upper  mantle  structure  in  the  same  regions,  bA  in  the  lower  mantle  part  of  the 

g 

ray  path,  and  b by  the  source  radiation  pattern.  All  these  factors  clearly  may  depend  upon  the 
source-receiver  configuration  through  source  take-off  angle  and  receiver  arrival  angle.  The 

Gutenberg-Richter  distance-depth  correction  may  be  considered  to  approximate  the  effects 
cs  us  IV1  ur  cr 

(b  + b + b + b + b ) globally;  the  biases  measured  then  really  measure  deviations  from 

g 

this  average  behavior.  The  term  b accounts  for  deviations  from  the  average  amplitude  in  a 

small  region  surrounding  the  source  due  to  radiation  pattern  effects.  Station  site  effects,  e.g., 

seismometer-ground  coupling  variations  due  to  the  medium  (hard  rock,  alluvium)  upon  which 

cr 

the  station  is  situated,  are  included  in  b 

Seismic  sources,  by  their  very  nature,  have  a tendency  to  be  located  in  regions  of  high 
lateral  inhomogeneity  and,  thus,  it  may  be  expected  that  features  such  as  the  anomalously  high 
attenuation  on  the  concave  side  of  island  arcs^  can  seriously  affect  determinations.  Studies 
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leading  to  the  development  of  plate  tectonic  theory  have  indicated  consistency  of  fault  plane  solu- 

s 

tions  and,  thus,  presumably  radiation  patterns  and  b over  large  source  regions. 

The  sources  used  in  Sec.  A have  been  separated  into  11  major  source  regions,  shown  in 
Fig.  1-4,  and  biases  computed  in  the  same  manner  as  before  for  events  from  each  of  these  source 
areas  reported  at  a given  station.  A value  of  mean  bias  is  not  calculated  unless  there  are  at 
least  100  events  in  any  particular  source  region  data  sample. 

For  16  of  the  72  stations  of  Sec.  A,  the  mean  bias  for  any  one  particular  source  region  dif- 
fers by  more  than  0.2  units  from  the  mean  bias  for  all  source  regions  given  in  Table  1-1. 

For  three  stations  in  highly  heterogeneous  seismic  regions  (PMG,  New  Guinea;  PMR,  Alaska; 
and  TSK,  Japan),  it  differs  from  the  mean  bias  for  all  source  regions  by  more  than  0.3  m^ 
units.  It  must  be  noted  here  that  for  some  stations  there  are  an  insufficient  number  of  events 
for  individual  source  region  biases  to  be  calculated  and,  thus,  comparison  of  biases  from  dif- 
ferent regions  is  severely  limited  to  the  better  stations.  For  some  of  the  best  stations,  e.g., 
TFO,  UBO,  W MO  which  have  reported  many  events  from  various  source  regions,  the  variation 
of  bias  with  source  region  is,  however,  less  than  0.2  units. 

Many  of  the  stations  considered  here  happen  to  be  concentrated  in  regions  which  are  small 
in  extent  compared  to  source-receiver  distances  and  we  may  therefore  expect  to  see  consisten- 
cies in  the  variations  of  mean  bias  with  source  region.  Three  particularly  small  receiver  re- 
gions are  those  containing  stations  in  Germany,  East  Africa,  and  the  western  U.S.  Table  1-3 
gives  the  variations  of  bias  with  source  region  for  each  of  these  areas. 

It  can  be  seen  that,  even  over  these  small  receiver  regions,  there  is  remarkably  little 
consistency  in  variations  of  mean  bias  with  source  region,  the  only  possible  exception  being 
that  for  events  in  region  4 (Japan)  to  the  western  U.S.  Similar  tables  for  other  receiver  areas 
such  as  Scandinavia  and  India  also  fail  to  reveal  any  correlation  of  bias  with  source  region. 

Had  we  selected  smaller  source  regions,  it  is  likely  both  that  regional  variations  in  bias  and 
their  consistency  across  receiver  areas  would  have  been  more  pronounced,  particularly  in  the 
case  of  ray  paths  traveling  down  descending  lithospheric  slabs.  Unfortunately,  the  data  base 
is  insufficient  to  test  this. 

c s us  M 

This  lack  of  correlation  may  be  taken  to  indicate  that  the  terms  (b  + b + b ),  however 
large  they  may  be  for  individual  ray  paths,  tend  to  average  out  in  such  a manner  that  they  can- 
not be  resolved  in  the  present  study.  It  also  indicates  that  the  Gutenberg-Richter  distance- 

depth  correlation  is  not  grossly  in  error.  In  particular,  much  of  the  observed  station  bias  may 
ur  c r 

be  due  to  the  term  b , and  its  variation  with  source  region  to  b . This  does  not  of  course 
cs  us  cr  ur 

imply  that  b and  b are  not  as  large  as  b and  b ; in  fact,  the  greater  lateral  heterogeneity 
in  source  regions  probably  means  that  they  will  be  larger,  but  they  cannot  be  resolved  by  the 
present  means.  R.  G.  North 

C.  TEMPORAL  VARIATIONS  IN  THE  SHAPE 
OF  THE  MAGNITUDE -FREQUENCY  CURVE 

Variations  in  the  shape  of  the  magnitude-frequency  curve  and,  in  particular,  the  slope  b of 
the  (dubiously)  linear  portion  of  the  curve,  where  b is  given  by  logN  = a — bm,  N being  the 
number  of  events  of  magnitude  (a  and  b constants),  have  been  interpreted^  as  indicative  of 
changes  in  the  state  of  stress  in  the  earth’s  crust  with  source  region  and/or  time.  In  a previous 
SATS,^  it  has  been  shown  that  these  variations  in  b are  of  dubious  statistical  significance. 


5 


TABLE  1-3 

VARIATIONS  OF  STATION  BIAS  WITH  SOURCE  REGION 
(With  Respect  to  the  Bias  Values  Listed  in  Table  1-1) 

Source  Region  (See  Fig.  1-4) 

- 

+ 0.01 

-0. 16 
-0.  02 
+ 0.  23 
+ 0. 08 

-0.09 

o 

-0.08 
-0.21 
+ 0.  07 
0.0 
-0.  07 

+ 0.01 
-0.03 
-0.21 
+ 0.05 

o 

*n  rx  o 

— o — 

o o o 

+ + + 

CO  — CO  00 

— 0 — 00 

o*  o*  o*  o'  o* 

+ 1 + + 

00 

o — 

— o 

o*  o* 

1 1 

rv 

»o 

CN 

O 

+ 

+ 0.  07 
0.0 

+ 0.05 

o 

+ 0.  02 
-0.04 
-0.  15 

-0.03 
-0. 13 

-0.  07 
-0.  07 
-0.06 
+ 0.  02 
+ 0.03 

in 

+ 0.01 
+ 0.18 
-0.09 
0.0 
+ 0.09 
+ 0.03 

-0.  07 
+ 0.05 
-0.  13 
-0. 25 
-0.03 

-0.09 
-0.01 
+ 0.04 
+ 0.09 
-0.04 
+ 0.  12 

o 

o 

+ 

-0.08 

-0.08 

-0.09 

-0.  13 

co 

-0.04 
+ 0.21 

+ 0.05 

+ 0.05 
+ 0.04 
-0.  15 
-0.06 

<N 

- 

o 

o* 

1 

Station 

Germany 

BNS 

CLL 

FUR 

GRF 

MOX 

STU 

East  Africa 
BHA 
CIR 
CLK 
KRR 

Western  U.  S. 
DUG 
EUR 
TFO 
TUC 
UBO 
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We  consider  here  the  effects  of  temporal  changes  in  the  set  of  stations  reporting  values 
in  the  ISC  Bulletins  for  1964-1973  upon  the  shape  of  magnitude -frequency  curves.  Such  changes 
may  be  expected  if  the  proportion  of  stations  reporting  which  have  large  negative  or  positive 
magnitude  biases  (see  Sec.  A,  above)  changes  with  time.  An  inspection  of  the  set  of  stations 
reporting  values  for  each  year  of  1964-1973  reveals  that  this  indeed  may  be  the  case.  For 
the  first  six  years  of  this  period,  a very  large  proportion  of  all  m^  reports  was  from  stations 
in  the  western  U.S.;  this  proportion  declined  considerably  over  1970-1973.  Since  low  upper 
mantle  Q beneath  the  western  U.S.  causes  substantial  negative  magnitude  biases  for  stations 
in  this  region  (see  previous  section),  this  implies  that  average  values  for  the  earlier  time 
period  will  be  lower  than  later  ones,  relative  to  "true”  m^.  Figure  1-5  shows  the  total  number 
of  station  reports  per  year  for  (a)  all  stations,  and  (b)  stations  in  the  western  U.S.  as  a func- 
tion of  year.  It  can  be  seen  that  the  proportion  of  station  m^  reports  from  western  U.S.  stations 
has  declined  from  over  60  percent  in  1964-1965  to  30  percent  in  1970-1973.  The  largest  change 
occurs  during  1969-1970. 

The  seismicity  studied  has  been  separated  into  that  for  1964-1969  and  1970-1973,  and  the 
annual  number  of  events,  in  m^  increments  of  0.1,  for  these  intervals  compared  with  that  for 
1964-1973.  Figure  I-6(a)  shows  ratios  of  annual  seismicity  for  each  of  the  time  periods  1964- 
1969  and  1970-1973  to  that  for  1964-1973.  It  is  clear  that  during  the  earlier  time  period  there 
were  many  more  events  reported  of  lower  magnitudes,  and  slightly  less  at  higher  magnitudes, 
than  in  the  later  time  period.  Application  of  the  station  biases  determined  in  the  previous  sec- 
tion as  corrections  to  station  magnitudes  before  calculation  of  average  [Fig.  1-6 ( b ) ] reduces 
the  difference  in  numbers  of  events  of  m^  > 4.6  between  the  two  time  periods  by  a factor  of  2. 
There  are  still  many  more  events  at  lower  magnitudes  for  the  earlier  time  period;  this  is  due 
to  the  remarkable  detection  performance  of  the  5 Vela  short-period  arrays  (BMO,  CFO,  TFO, 
UBO,  and  WMO).  These  stations  closed  down  (TFO,  WMO),  or  reported  a smaller  proportion 
of  later  events  during  1969-1970.  Obviously,  the  detection  capability  of  the  global  network  of 
stations,  as  measured  in  terms  of  m^,  has  been  seriously  degraded  by  the  removal  of  these 
Vela  stations.  It  is  surprising  that  although  there  are  more  stations  reporting  in  later  years, 
they  do  not  do  as  well  at  lower  magnitudes  as  the  fewer  stations  of  1964-1969. 

It  can  be  concluded  that  temporal  changes  in  station  distribution  can  cause  variations  in 
the  distribution  of  event  magnitudes  reported  by  a global  network.  Similar  regional  changes 
also  can  be  expected  as  the  location  of  stations  relative  to  world  seismicity  varies. 

R.  G.  North 

D.  COMPARISON  OF  RECENT  ESTIMATES  OF  MAGNITUDE  BIAS 

In  Sec.  A,  we  have  presented  a set  of  station  biases  determined  directly  from  station  ampli- 
tude readings.  It  is  interesting  to  compare  these  results  with  two  recent  studies  of  magnitude 
7 8 

bias  by  Marshall  and  Evernden.  Table  1-4  lists  the  biases  from  all  three  studies  for  a set  of 
stations  common  to  at  least  two  studies.  Results  from  Marshall  and  Evernden  have  each  been 
subjected  to  a base  line  correction  so  that  station  COL  (College,  Alaska)  has  the  value  +0.01  in 
each  case. 

7 

The  results  of  Marshall  (kindly  supplied  to  us  in  advance  of  publication)  were  obtained  using 
estimates  of  the  mean  attenuation  coefficient  Q of  the  crust  and  upper  mantle  beneath  the  re- 
ceiver obtained  from  observations  of  velocity.  The  general  agreement  with  the  Lincoln 


7 


TABLE  1-4 

COMPARISON  OF  ESTIMATES  OF  STATION  MAGNITUDE  BIAS 

Station 

Magnitude  Biases 

Station 

Magnitude  Biases 

Table 

1-1 

Marshall  ^ 

Evernden^ 

Table  1-1 

Marshal  1^ 

Evernden® 

TFO 

-0. 

32 

-0.21 

-0. 25 

COLt 

+ 0.01 

+ 0.01 

+ 0.01 

BMO 

-0. 

29 

-0.21 

-0. 25 

KEV 

+ 0.  02 

+ 0.06 

+ 0.37 

GOL 

-0. 

28 

-0. 23 

LPS 

+ 0.04 

+ 0.04 

ALQ 

-0. 

20 

-0. 23 

-0.  17 

HFS 

+ 0.05 

+ 0.03 

WMO 

-0. 

17 

-0.18 

KOD 

+ 0.06 

-0.  12 

DUG 

-0. 

15 

-0.  28 

-0.31 

TRN 

+ 0.07 

+ 0.  15 

NOR 

-0. 

14 

-0.12 

+ 0.20 

KON 

+ 0. 07 

-0.05 

+ 0.  16 

UBO 

-0. 

11 

-0.21 

+ 0.03 

KHC 

+ 0.  10 

+ 0. 10 

LAO 

-0. 

10 

+ 0.01 

PMG 

+ 0.10 

+ 0.21 

WIN 

-0. 

09 

+ 0.  12 

SHL 

+ 0.11 

-0.  02 

NAO 

-0. 

09 

-0.05 

KJN 

+ 0.14 

+ 0.06 

BUL 

-0. 

07 

-0.01 

POO 

+ 0.  17 

+ 0.  15 

+ 0. 22 

PRE 

-0. 

07 

+ 0.23 

HYB 

+ 0. 19 

+ 0.  15 

CPO 

-0. 

07 

+ 0. 10 

+ 0.06 

SJG 

+ 0.  24 

+ 0.  24 

BOZ 

-0. 

06 

-0.21 

STU 

+ 0. 29 

+ 0.01 

+ 0.27 

ALE 

-0. 

04 

-0.  12 

NDI 

+ 0. 33 

+ 0.  15 

+ 0.  18 

EKA 

0 

+ 0.01 

EDM 

+ 0.37 

+ 0.  10 

t Reference  station. 
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Laboratory  results  is  excellent  (differences  in  estimated  bias  seldom  exceed  0.1  m^  units).  It 
is  interesting  to  note,  however,  that  the  total  range  of  biases  given  by  Marshall  is  approximately 
0.4  units,  compared  with  a range  of  nearly  0.7  units  in  the  present  study. 

Marshall's  results  raise  the  possibility  that  there  may  be  some  correlation  between  station 
magnitude  bias  and  local  Pn  velocity.  To  test  this.  Fig.  1-7  shows  a plot  of  bias  against  Pn 
velocity  for  40  stations.  The  Pn  values  are  those  obtained  by  Marshall  (private  communication). 
A definite  trend  is  seen,  which  may  be  approximated  by  the  expression 

bias  = (P  velocity)  - 8.05 

Because  of  the  scatter  in  the  data  points,  we  do  not  suggest  that  this  relation  has  a useful  pre- 
dictive capability.  It  is  interesting,  however,  and  gives  some  support  to  the  arguments  of 

7 

Marshall. 

g 

In  a second  recent  study,  Evernden  has  estimated  station  biases  using  a relationship  be- 

9 

tween  bias  and  background  noise  level  presented  in  an  earlier  paper.  The  agreement  of  these 
results  with  the  first  two  studies  is  reasonably  good,  in  most  cases. 

In  addition  to  Pn  velocity  and  background  noise  level,  a third  parameter  which  might  be 
expected  to  correlate  with  station  magnitude  bias  is  the  station  travel  time  correction.  Fig- 
ure 1-8  shows  a plot  of  magnitude  bias  against  the  station  corrections  listed  by  Herrin  and 
1 0 

Taggart.  Quite  clearly,  there  is  no  correlation  between  these  parameters. 

M.  A.  Chinnery 
R.G.  North 


E.  SEISMIC  MAGNITUDE  SIMULATION  STUDIES 

A statistical  model  for  the  seismic  magnitudes  measured  at  a network  of  stations  was 
formulated  in  our  last  SATS.  Some  implications  of  the  model  were  developed  as  were  algo- 
rithms for  the  estimation  of  model  parameters  and  magnitude.  The  basic  model  continues  to 
be  used  to  interpret  magnitude  data, and  it  was  felt  that  computer  simulations  should  be  done  to 
allow  various  controlled  experiments  in  which  all  model  parameters  and  event  magnitudes  were 
known.  Such  a simulation  may  indicate  how  the  model  might  be  reformulated  to  better  represent 
actual  data,  and  we  plan  to  apply  proposed  estimation  algorithms  to  simulated  data  to  evaluate 
their  performance. 

A simulation  program  has  been  written  which  is  in  terms  of  individual  station  magnitudes 
and  does  not  take  account  of  geographical  distribution  of  either  stations  or  events.  The  simula- 
tion generates  a sequence  of  true  magnitude  values  and  for  each  of  these  randomly  generates  the 
magnitudes  seen  at  a network  of  stations.  The  user  specifies  the  number  of  stations  which  must 
detect  to  declare  an  event  detected  by  the  network.  The  program  outputs  the  simulated  magni- 
tude data  as  well  as  the  true  magnitude  for  each  event  detected.  All  this  is  driven  by  a pseudo- 
random number  generator  subroutine.  The  main  input  options  are  as  follows. 

1.  True  Magnitude  Distributions 

(a)  The  user  can  select  to  generate  true  magnitude  values  by  sampling  from 
a probability  density  function  which  is  of  the  form  Ce  on  an  interval 
from  m = mQ  to  m = 7.0.  C is  automatically  selected  to  properly  nor- 
malize the  distribution.  The  user  can  choose  the  minimum  magnitude 
mQ  and  the  slope  ft  of  the  seismicity  curve.  This  choice  for  the 
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generation  of  true  magnitude  corresponds  to  the  common  assumption  of 
a linear  relation  between  average  number  of  events  and  magnitude.  The 
nominal  assumption  of  a factor  of  10  change  in  the  number  of  events  for 
each  magnitude  unit  is  equivalent  to  using  /?  = In  10  =*  2.3. 

(b)  A second  alternative  is  that  all  events  have  the  same  true  magnitude.  In 
this  case,  the  user  selects  the  value. 

(c)  The  third  alternative  is  for  the  user  to  supply  a special  subroutine  which 
allows  him  to  specify  any  probability  density  function  at  all  for  the  true 
magnitudes. 

2.  Bias  and  Scatter  Parameters 

Given  true  magnitude  m the  predetection  magnitude  at  station  i is  given 
by: 

m . - m + B.  + e.  (I- 1 ) 

ill 

where  B.  is  a constant  bias  and  e.  is  a zero-mean  Gaussian  random  vari- 
l l 

able  with  standard  deviation  o^.  The  user  must  specify  the  B^  and  a. 
values. 


3.  Station  Detection  Characteristics 


(a)  The  basic  option  for  the  detection  characteristics  of  a station  is  that  the 
probability  of  detection  is: 

/ m . - G? 

Pr  [detection  /m^]  = $1  — - — - — — 

\ y* 

where  $ is  the  error  function,  G*  is  the  50-percent  detection  threshold, 
and  y*  determines  the  rate  at  which  the  probability  goes  from  zero  to 
one.  The  user  supplies  the  G*  and  y*.  For  a given  m^,  the  simulation 
randomly  (using  the  proper  probability)  decides  if  a detection  has  oc- 
curred. It  then  reports  either  m-  or  the  "not  seen"  value. 

(b)  In  recognition  that  the  error  function  may  not  be  a satisfactory  detection 
curve  in  all  cases,  the  program  allows  the  user  to  supply  a subroutine 
which  can  modify  the  basic  detection  curve.  The  problem  of  saturation 
of  stations  for  large  events  can  be  conveniently  accommodated  using  this 
option. 


(1-2) 


One  imminent  application  of  the  magnitude  simulation  program  is  to  check  for  the  possible 
influence  of  station  detection  characteristics  upon  the  station  bias  estimates  obtained  by  North 
(see  Sec.  A).  Given  North's  station  bias  estimates  plus  estimates  of  other  model  parameters, 
we  can  generate  simulated  magnitudes  for  a network.  The  algorithms  used  by  North  to  estimate 
biases  can  then  be  applied  to  these  data  to  see  if  the  now  known  biases  built  into  the  model  will 
be  regenerated  by  the  estimation  algorithms. 

The  model  parameters  required  to  complete  such  an  experiment  are  the  seismicity  param- 
eter /3,  the  scattering  parameter  o.,  and  the  detection  parameters  G*  and  y ? . An  effort  has 
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been  made  to  obtain  estimates  of  the  G?  and  y*  which  could  be  used  for  the  experiment.  These 
estimates  are  probably  not  optimal  but  should  be  adequate  for  the  application  at  hand. 

A log-likelihood  function  for  the  estimation  of  seismicity  and  detection  parameters  of  a 
single  station  in  terms  of  the  reported  amplitudes  and  periods  at  a station  was  given  in  our  last 
SATS.  In  terms  of  station  magnitudes,  this  becomes 


Ni 

lnLi=  i 

j=l 


In  p + P(G*  - | 07? 


2 


) + In 


(1-3) 


where  i indicates  the  station,  N.  is  the  number  of  events  detected  at  station  i,  and  j is  the  index 

1 J 

on  the  events.  We  have  written  and  used  a program  to  search  the  /3,  G? , y*  space  for  the  maxi- 
mum value  of 

c = 2>Li 

and  to  thus  estimate  0,  G*,  and  y*.  For  the  results  given  below,  the  grids  for  this  search 
were: 


(3  = 1.5  to  2.4  in  steps  of  0.1 


G^  = 3.4  to  5.8  in  steps  of  0.1 
y*  = 0.1  to  0.55  in  steps  of  0.05 


A preliminary  simulation  for  15  stations  and  7 500  detected  events  was  run  with  seismicity 
parameter  (3  = 2.3,  o\  = 0.3  for  all  stations,  and  other  parameters  as  given  in  Table  1-5.  Pa- 
rameters were  selected  arbitrarily.  The  object  was  to  check  the  estimation  algorithms  with 
simulated  data  using  known  parameter  values.  Given  the  quantizations  used  for  estimation  and 
the  small  number  of  events,  performance  was  quite  satisfactory.  The  seismicity  parameter 
was  underestimated  as  2.1  rather  than  the  correct  2.3  value.  Estimates  of  G?  and  y*  are  also 
good  with  some  tendency  to  underestimate.  The  bias  at  some  stations  caused  no  problem.  The 
quality  of  estimates  generally  decreased  with  the  number  of  magnitude  reports  available.  In 
anticipation  of  results  with  real  data,  we  note  that  when  we  assigned  (3  = 1.6  for  estimation  pur- 
poses (although  the  true  value  was  2.3),  the  resulting  detection  parameter  estimates  did  not 
change  very  much.  As  might  be  expected,  the  estimated  50-percent  detection  thresholds  de- 
creased when  (3  was  fixed  at  the  lower  value. 

Table  1-6  shows  the  results  obtained  using  real  data.  Based  upon  numbers  of  reports  to 
the  ISC,  North  selected  a set  of  72  stations  to  study  in  detail.  We  used  the  data  from  those  sta- 
tions from  1964  through  197  3.  No  data  from  1969  were  used  due  to  a technical  problem.  The 
number  of  available  reports  per  station  ranged  from  1106  for  WIN  to  27,915  for  BMO.  The 
search  over  (3  as  well  as  detection  parameters  gave  the  unexpected  low  value  of  (3  = 1.6.  The 
estimated  50-percent  detection  magnitudes  ranged  from  3.8  to  5.4  and  the  scatter  parameters 
cluster  near  0.3.  The  estimates  of  50-percent  threshold  increased  by  0.0  to  0.2  magnitude  unit 
when  /3  was  fixed  at  the  value  2.3. 

Some  additional  work  is  required  but  we  believe  that  we  now  have  estimates  of  most  param- 
eters needed  so  that,  given  North's  bias  estimates,  we  can  undertake  a meaningful  simulation. 
Our  next  step  will  be  to  verify  within  the  context  and  assumptions  of  our  model  that  the  station 
bias  estimates  obtained  by  North  have  not  been  excessively  influenced  by  station  detection 
characteristics.  R.  T.  Lacoss 

M.  A.  Chinnery 
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TABLE  1-5 

ESTIMATE  OF  DETECTION  CHARACTERISTICS 
AND  SEISMICITY  FROM  SIMULATED  DATA  FOR  15  STATIONSt 


Stotian 

Number 
af  Reports 

Model  Parameters 

Estimoted  Parameters 

P = 

2.  1 

p = 

1.6 

Bios 

G* 

i 

G* 

i 

r: 

G* 

i 

1 

5451 

-0.3 

3.8 

0.2 

3.8 

0.  15 

3.7 

0.  15 

2 

4350 

-0.3 

3.9 

0.2 

3.9 

0.20 

3.8 

0.  15 

3 

3510 

1 

© 

CO 

4.0 

0.2 

4.0 

0.20 

3.9 

0.  15 

4 

3582 

-0.2 

4.1 

0.2 

4.1 

0.20 

4.0 

0.  15 

5 

2886 

CN 

O 

1 

4.2 

0.2 

4.2 

0. 20 

4.1 

0.  15 

6 

2293 

-0.2 

4.3 

0.2 

4.3 

0. 20 

4.2 

0.  15 

7 

2221 

-0.1 

4.4 

0.2 

4.4 

0.  15 

4.3 

0.  10 

8 

1891 

-0.1 

4.5 

0.2 

4.5 

0.20 

4.4 

0.  15 

9 

1512 

-0.1 

4.6 

0.2 

4.5 

0.  15 

4.5 

0.  15 

10 

1513 

0.0 

4.7 

0.2 

4.7 

0. 20 

4.6 

0.  15 

11 

1228 

0.0 

4.8 

0.2 

4.8 

0.  15 

4.7 

0. 15 

12 

995 

0.0 

4.9 

0.2 

4.9 

0.  15 

4.8 

0. 15 

13 

743 

0.0 

5.0 

0.2 

4.9 

0. 10 

4.9 

0. 15 

14 

646 

0.0 

5.  1 

0.2 

5.0 

0.  10 

5.0 

0.  10 

15 

621 

0.1 

5.2 

0.2 

5.1 

0. 10 

5.1 

0.  10 

t Model  parameters  were  o\  = 0.3,  P = 2. 3 and  Bios,  G^,  yjf  os  shown  in  the  toble. 
Estimoted  detection  parameters  G-*,  y*'  shown  far  estimoted  seismicity  parameter 
P = 2.  1 and  far  P improperly  constrained  to  p = 1.6  for  estimotion  purposes. 
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ESTIMATES  OF 

TABLE  1-6 

DETECTION  CHARACTERISTICS  FOR  72  STATIONSt 

Station 

P = 

1.6 

P = 

2.3 

Stotion 

p = 

1.6 

-p  = 

2.3 

G* 

i 

G* 

i 

?? 

G? 

i 

G.* 

i 

r; 

CPO 

4.2 

0.30 

4.3 

0.  35 

KJF 

4.8 

0. 25 

4.9 

0.30 

TUL 

4.7 

0. 35 

4.8 

0.35 

KJN 

4.7 

0.  25 

4.8 

0.30 

WMO 

4.  1 

0.25 

4. 1 

0. 25 

KON 

5.0 

0. 30 

5.0 

0. 25 

ALQ 

4.4 

0. 35 

4.5 

0. 35 

NAO 

3.9 

0.35 

4.0 

0. 35 

BKS 

5.4 

0.50 

5.6 

0. 55 

NUR 

4.8 

0.30 

4.9 

0.30 

BMO 

3.9 

0. 25 

4.0 

0.  25 

EKA 

4.6 

0.30 

4.7 

0.30 

BOZ 

4.5 

0.30 

4.6 

0. 30 

KHC 

4.8 

0. 25 

4.9 

0. 25 

DUG 

4. 1 

0.30 

4.2 

0.30 

KRA 

5.  1 

0. 25 

5.2 

0. 25 

EUR 

4.2 

0.30 

4.3 

0. 30 

LJU 

5.2 

0.30 

5.4 

0. 35 

GOL 

4.3 

0.30 

4.4 

0. 30 

LOR 

4.8 

0.50 

4.8 

0.50 

LAO 

4.0 

0.30 

4.  1 

0.30 

NIE 

4.7 

0.20 

4.8 

0.20 

LON 

4.6 

0.30 

4.7 

0.30 

PRU 

5.0 

0.30 

5.0 

0.30 

NEW 

4.8 

0.35 

4.9 

0. 35 

BNS 

5.2 

0.30 

5.3 

0.30 

TFO 

3.8 

0. 25 

3.9 

0. 30 

CLL 

4.9 

0.  25 

5.0 

0. 25 

TUC 

4.6 

0.30 

4.7 

0. 30 

FUR 

5.0 

0.35 

5.  1 

0. 35 

UBO 

3.9 

0. 35 

4. 1 

0.30 

GRF 

5.  1 

0.  25 

5.2 

0.25 

BHA 

4.4 

0. 25 

4.5 

0.  25 

MOX 

4.8 

0.30 

4.9 

0.30 

BNG 

4.2 

0. 25 

4.3 

0. 25 

STU 

5.2 

0. 30 

5.3 

0.30 

BUL 

4.5 

0.20 

4.6 

0. 25 

COL 

4.5 

0.35  ‘ 

4.6 

0.30 

CIR 

4.4 

0.25 

4.6 

0.30 

GIL 

4.5 

0.35 

4.6 

0.35 

CLK 

4.4 

0.25 

4.5 

0.25 

PMR 

4.6 

0.30 

4.7 

0. 35 

KRR 

4.4 

0.20 

4.5 

0. 20 

GDH 

5.0 

0. 35 

5. 1 

0. 35 

PRE 

4.6 

0.25 

4.8 

0.30 

KTG 

4.9 

0.35 

5.0 

0.  35 

WIN 

4.6 

0. 25 

4.7 

0. 25 

NOR 

4.6 

0.25 

4.7 

0.  25 

HYB 

5.  1 

0.35 

5.2 

0.30 

ALE 

4.9 

0.30 

5.0 

0.30 

KOD 

4.9 

0.30 

5. 1 

0. 35 

MBC 

4.8 

0. 40 

4.9 

0. 40 

NDI 

5.2 

0.35 

5.3 

0. 30 

NP- 

4.3 

0. 30 

4.4 

0.30 

POO 

5.0 

0.50 

5.2 

0. 55 

RES 

5.0 

0.40 

5.1 

0.35 

SHL 

4.7 

0.30 

4.8 

0. 30 

EDM 

5.3 

0. 35 

5.3 

0.35 

ASP 

4.6 

0.30 

4.7 

0.  30 

PNT 

4.9 

0.25 

5.0 

0.  25 

CAN 

4.7 

0.30 

4.8 

0.30 

TSK 

4.7 

0.40 

4.9 

0. 45 

MUN 

5.0 

0.30 

5.1 

0.30 

TRN 

4.9 

0.  25 

5.0 

0.30 

RIV 

5.3 

0.30 

5.4 

0. 30 

LPS 

4.8 

0.30 

4.9 

0.30 

COP 

5.4 

0. 25 

5.5 

0. 25 

SJG 

4.9 

0.35 

5.0 

0.  30 

HFS 

4.2 

0. 35 

4.3 

0. 35 

CAR 

4.9 

0.25 

5.0 

0.25 

KEV 

4.9 

0.25 

5.0 

0.  25 

PMG 

4.8 

0.25 

4.9 

0. 25 

t Computed  using  station  magnitude  reports  from  1964  through  1973  except  for  1969.  Estimotes 
of  50-percent  threshold  G?  ond  spread  foctors  yj  shown  for  seismicity  parameter  p = 1.6  which 
resulted  from  searching  over  p ond  for  p constrained  to  2.3  for  estimotion  purposes. 
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F.  MAGNITUDE  DIFFERENCES  BETWEEN  STATION  PAIRS 

In  the  previous  work,  an  expression  was  given  for  the  expected  value  of  the  magnitude 
observed  at  station  2,  given  the  magnitude  m^  observed  at  station  1.  This  may  be  conveniently 
rewritten  in  terms  of  the  magnitude  difference  m^  — m^,  as  follows: 


2 , 2 

2 ® 

E[(m2  -m1)/m1]  = B2  - -/3ct1  + 


+a 


2 + >2 


X Z 


ml  + B2~B1  ~GZ 


I 2” 2 

J G 1 + a2  + Y 


*2 

2 


(1-4) 


This  formula  is  derived  assuming  a linear  frequency  magnitude  relation  of  the  form 

In  N = a — /3m  . (1-5) 

1 

The  parameters  involved  are  as  follows: 

B.  = bias  of  station  i 

1 

2 

(J ^ = variance  of  magnitude  readings  at  station  i 

G*  = detection  threshold  of  station  i 
2 

y . = variance  of  detection  probability  curve  for  station  i. 

1 

The  function  z(x)  is  defined  as  previously. 

Equation  (1-4)  is  interesting  for  several  reasons.  It  shows  that  the  magnitude  difference 
(m^  — m^),  averaged  at  each  value  of  m^,  has  an  expected  value  independent  of  the  detection 
characteristics  of  m^.  Furthermore,  the  appearance  of  the  quantity  p in  the  expression  sug- 
gests that  information  about  the  seismicity  distribution,  Eq.  (1-5),  can  be  obtained  even  though 
the  data  used  in  the  evaluation  of  Effm^  — m^)/m^J  are  not  complete. 

As  a test  of  the  usefulness  of  Eq.  (1-4),  we  report  here  the  results  of  plotting  the  average 
value  of  (m^  — m^)  against  for  several  station  pairs.  The  data  used  were  those  of  the 
Bulletin  of  the  International  Seismological  Center  (ISC),  for  the  years  1966-1970,  which  are 
available  on  tape.  For  a given  station  pair,  events  were  selected  that  had  a recorded  magnitude 
at  each  station,  and  which  lay  in  the  distance  range  30°  to  90°  from  each  station  (to  minimize 
possible  errors  in  the  distance  correction). 

An  example  of  the  resulting  data  is  shown  in  Fig.  1-9,  for  the  station  pair  TFO  — UBO,  using 
reference  station  UBO.  The  points  show  the  average  value  of  the  magnitude  difference 
mTFO  “ mUBO'  at  increments  0*1  m]3  units  in  m^rgQ*  The  vertical  bars  indicate  2 standard 
deviations  on  either  side  of  the  mean  (i.e.,  9 5-percent  confidence  bounds).  Approximately 
600  measurements  were  used  in  the  determination  of  each  point  in  the  range  4.0  < mUBO^ 

Fewer  observations  were  available  outside  this  range,  and  this  is  reflected  in  the  error  bounds 
associated  with  each  mean. 

The  solid  curve  shows  the  prediction  of  Eq.  (1-4)  using  the  following  parameters: 

B,  - B,  - Ba2  = -0.40 
2 1 1 

+ °2  ~ 0**9 
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G2  = 3.68 
y2  - °-32  • 

A very  good  fit  can  be  obtained,  though  some  questions  about  the  uniqueness  of  these  parameters 
remain. 

A similar  curve  can  be  obtained  from  the  same  two  stations,  but  using  TFO  as  a reference. 
This  is  shown  in  Fig.  I- 10.  The  curve  through  the  data  corresponds  to  the  parameters: 

B,  - B,  - Ba,2  = -0.01 

Z 1 1 

CT2  + <72  = 0.19 

1 z 

G2  = 4.15 

y £ - 0.50 

The  bias  levels  in  Figs.  1-9  and  -10  are  consistent  with  the  following: 

bubo  ” btfo  = 0-21 

aUBO  = °-30 

aTFO  = 0,31 
P = 2.30 

The  difference  in  biases  of  0.21  is  in  excellent  agreement  with  the  difference  in  bias  found  by  a 
quite  different  method  elsewhere  in  this  report  (Sec.  A). 

Provisional  conclusions  from  these  results  are  that  the  theoretical  model  appears  to  fit 
the  observed  data  well,  in  certain  cases.  The  assumption  of  linear  seismicity  appears  to  be 
valid.  Note  that  the  value  /3  = 2.30  corresponds  to  a "b-value"  of  1.0.  In  cases  like  those 
shown,  it  appears  possible  to  derive  valuable  information  on  detection  parameters  and  biases. 

Other  station  pairs  do  not  always  behave  so  well.  Two  extreme  examples  are  shown  in 
Figs.  I- 11  and  -12  for  the  station  pairs  AIjQ  — UBO  and  EUR  — UBO.  These  curves  do  not  show 
a well-defined  base  level.  This  is  interpreted  as  being  due  to  the  effects  of  instrument  satura- 
tion, which  are  overlapping  with  the  effects  of  detection  level  in  the  range  5.0  ^ ^ 5.5.  Ef- 

forts to  extract  useful  parameters  in  cases  of  this  type  are  continuing. 

M.  A.  Chinnery 
R.  T.  Lacoss 


G.  ESTIMATION  OF  SEISMIC  SOURCE  AT  1 Hz  FROM  log(A/T) 

Short-period  magnitudes  have  traditionally  been  calculated  by  = log(A/T)  + Q(A,  h), 
where  A is  the  waveform  amplitude,  peak  to  trough,  corrected  for  instrument  response  at 
period  T,  and  Q(A,  h)  is  the  Gutenberg-Richter  correction  for  geometric  spreading  as  a function 
of  epicentral  distance  A and  depth  of  focus  h. 

Log  (A/T)  after  instrument  correction  is  usually  interpreted  to  be  a measure  of  the  source 
energy  at  1 Hz  arriving  at  the  seismometer.  This  is  approximately  true, but  the  effect  of  source 
function  and  instrument  response  on  the  time-domain  waveform  can  produce  a bias  which  varies 
with  the  corner  frequency  of  the  source  model  and  the  type  of  high-frequency  decay,  e.g.,  u>- 
squareor  a? -cube.  Furthermore,  earth  attenuation  which  is  frequency  dependent  has  an  effect 
which  is  not  properly  incorporated  in  the  Q(A,  h)  term  for  obtaining  m^  from  log  (A/T). 
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Figures  1-13,  -14,  and  -15  show  curves  of  corrected  log  (A/T)  for  three  source  models, 

11  1213 

recorded  on  the  standard  short-period  Vela  seismometer.  The  Sharpe  and  Brune  models  ' 

-2  14  -3 

decay  as  w at  high  frequencies,  and  the  Savage  model  decays  as  cj  . The  spectra  for 

these  models  are  scaled  so  that  at  0 Hz,  the  amplitude  equals  fQ^,  where  f is  the  corner  fre- 
quency. The  spectra  and  time-domain  waveforms  are  shown  in  Figs.  1-21,  -22,  and  -23. 

In  each  of  Figs.  1-13,  -14,  and  -15,  the  displacement  source  at  1 Hz  is  shown  as  a dashed 

line.  The  three  curves  show  the  effect  of  including  earth  attenuation,  using  a causal  attenuation 

15 

operator  given  by  Carpenter.  From  these  figures  it  is  clear  that  there  is  a significant  bias 
of  log  (A/T)  about  the  source  displacement  at  1-sec  period,  which  varies  with  attenuation  and 
corner  frequency  of  the  source  model.  Curiously,  the  attenuated  source  with  a t*  = 0.2  5,  when 
passed  through  the  seismometer  produces  an  A/T  value,  which  when  corrected  for  Vela  seis- 
mometer response  is  closest  to  the  unattenuated  source  spectrum  at  1 Hz,  over  a wide  range  of 
corner  frequencies. 

* The  bias  of  corrected  values  of  log  (A/T)  from  the  source  spectrum  at  1 Hz  is  displayed  in 

Fig.  1-16  for  each  source  model  and  several  values  of  t * as  a function  of  corner  frequency. 

-2  -3 

One  sees  that  the  biases  for  the  u:  and  u;  models  are  rather  similar  for  corner  frequencies 

down  to  ~0.25  Hz. 

One  can  calculate  mfe  for  corner  frequencies  f of  explosions  by  using  empirical  relations 
for  yield  Y in  kilotons,  equivalent  cavity  radius  a in  meters  and  m^  (Ref.  16).  If  we  assume  that 

a ~ 100  Y1/3 
mb  ~ 3.8  + log10  Y 


and 


f = -2— 

° 

in  the  Sharpe  model,  then,  we  cap  calculate  as  a function  of  a.  This  gives 
mb  ~ 3 log  (167/fQ) 

where  we  have  assumed  that  a = 5 km/sec.  Values  of  explosion  are  tabulated  between  the 

biases  for  the  Sharpe  and  Brune  models  in  Fig.  1-16.  These  m^  values  are  rough  averages,  so 

that  individual  events  may  not  have  corner  frequencies  consistent  with  the  m^  scale  shown. 

Also  shown  in  Fig.  1-16  are  the  estimates  of  corner  frequencies  for  the  explosion  Milrow  and 

four  Aleutian  earthquakes  measured  by  Wyss,  Hanks,  and  Liebermann  from  long-period  tele- 
17 

seismic  records.  It  is  apparent  that  the  corner  frequencies  for  these  earthquakes  are  an 
order  of  magnitude  lower  than  for  Milrow  for  roughly  comparable  m^. 

For  smaller  magnitude  explosions  between  4 and  5,  it  appears  that  the  bias  may  be 
higher  than  for  earthquakes  of  similar  magnitude,  along  similar  transmission  paths.  This 
should  produce  a measurable  effect  on  the  M discrimination  criterion  at  low  magnitudes. 

The  bias  described  above  should  vary  with  instrument  response.  We  plan  to  calucate  sim- 
ilar biases  for  WWSSN  and  SRO  recording  systems,  and  determine  the  necessary  corrections 
for  log  (A/T)  to  obtain  consistent  estimates  of  the  source  spectrum  at  1 Hz. 

C.  W.  Frasier 
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H.  EVIDENCE  FOR  a; -CUBE  SCALING  FROM  SHORT-PERIOD  AMPLITUDES 

AND  PERIODS  OF  THE  RAT  ISLAND  SEQUENCE  (196  5) 

We  are  studying  the  effect  of  seismic  scaling  directly  in  the  time  domain  by  analyzing  the 
amplitudes  and  periods  of  teleseismic  P-waves  from  moderate-sized  earthquakes  of  the  Rat 

Island  sequence.  The  motivation  for  this  study  is  the  clear  linear  trend  of  log  amplitude  with 

1 

period  for  this  sequence  reported  by  USGS  EDR.  In  the  previous  SATS,  we  showed  the  trend 
of  log  (A/T)  vs  T for  over  300  events  recorded  at  four  Vela  arrays.  Since  then  we  have  ex- 
tended the  data  base  to  868  events  recorded  at  five  Vela  arrays  and  four  WWSSN  stations.  We 

have  also  generated  synthetic  records  of  the  far-field  displacement  including  instrument  re- 

11  12  13 

sponse  and  earth  attenuation  for  two  lo  -square  models  proposed  by  Sharpe  and  Brune, 

14 

and  an  a? -cube  model  by  Savage. 

Under  the  assumption  of  geometric  similarity  of  small-  and  moderate-sized  earthquakes, 
it  can  be  seen  that  predicted  variations  of  log  (A/T)  with  T are  quite  different  for  a? -square  and 
a; -cube  models  when  the  corner  frequency  moves  through  the  short-period  band.  Using  the 
Vela  station  data,  we  can  demonstrate  that  the  a? -square  models  may  be  rejected  in  favor  of 
the  w-cube  model  suggested  by  Savage. 

The  recording  stations  and  the  aftershock  area  are  shown  in  Fig.  1-17.  The  Vela  arrays 
arc  BMO,  UBO,  TFO,  WMO,  and  CPO.  They  fall  in  the  azimuthal  sector  150°  to  170°  from 
the  source  region,  and  within  the  distance  range  44°  to  68°.  WWSSN  stations  EUR  and  ALQ 
are  located  within  the  above  sector,  and  stations  COL  and  NUR  outside  the  sector. 

Figure  I- 18(a)  shows  the  scatter  of  reported  amplitudes  with  periods  at  the  Vela  station 
BMO  with  a histogram  of  the  amplitudes  in  Fig.  l-18(c).  In  Fig.  1 - 1 8 ( b ) the  mean  and  standard 
deviation  of  the  log  amplitude  values  are  displayed  as  a function  of  T.  The  linear  trend  of 
log  A with  period  is  quite  striking  and  is  typical  of  the  data  from  the  other  Vela  stations.  In 
Fig.  I-18(d),  log  (A/T)  vs  period  is  displayed,  and  the  clear  trend  with  T is  still  apparent 
although  the  slope  is  smaller.  In  the  figures  which  follow,  we  plot  log  (A/T)  rather  than  log  A 
simply  because  it  is  a measure  of  m^,  and  we  shall  attempt  to  relate  m^  to  corner  frequency 
for  these  data. 

The  log  (A/T)  values  for  the  four  other  Vela  stations  are  shown  in  Fig.  1-19.  The  numbers 
for  each  station  equal  the  number  of  events  reported  during  the  time  period  scanned  in  the  EDR. 
Although  some  data  scatter  is  apparent,  the  mean  values  of  log  (A/T)  clearly  increase  with  T 
over  a wide  range  of  periods  with  similar  slopes  at  all  Vela  stations. 

Figure  1-20  contains  similar  plots  of  amplitudes  reported  by  four  WWSSN  stations.  In  this 
figure,  we  see  an  increase  in  log  (A/T)  with  T at  EUR  and  NUR,  but  no  reliable  increase  at 
ALQ  and  COL.  Comparing  this  figure  to  Fig.  1-19,  one  can  see  several  important  differences 
between  reported  amplitudes  at  Vela  and  WWSSN  stations.  First,  the  number  of  events  reported 
by  Vela  stations  is  2 to  10  times  the  number  reported  by  the  WWSSN  stations.  The  range  of  re- 
ported periods  is  wider  at  the  Vela  stations,  and  the  lower  limit  of  average  amplitudes  at  short 
periods  is  below  the  lower  limits  at  the  standard  stations.  In  terms  of  data  quality  and  consis- 
tency of  reported  amplitudes,  it  appears  that  the  Vela  arrays,  while  operating,  performed  much 
better  than  the  standard  stations.  This  is  discussed  in  another  section  on  station  biases  by  North. 

The  source  models  we  used  are  those  of  Sharpe/*  Brune/2,13  and  Savage/"*  We  are  pri- 
marily interested  in  scaling  effects  with  corner  frequency  f . Therefore,  we  ignore  all  fre- 
quency independent  factors  such  as  spherical  spreading  and  consider  only  the  far-field  displace- 
ment wave  shapes  and  Fourier  spectra,  including  phase.  The  variation  in  both  azimuth  and 
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take-off  angle  of  the  rays  as  they  leave  the  source  is  fairly  small,  except  toward  NUR  and  COL, 
so  we  shall  assume  that  the  variation  in  radiation  pattern  is  negligible  in  the  take-off  direction 
to  the  Vela  arrays. 

The  Fourier  transforms  of  the  scaled  source  functions  are,  in  chronological  order. 


Sharpe:  U^f)  = fQ'3/{  1 - (f/fQ)2  + (f/fQ)} 

(1-6) 

Brune:  U,(f)  = fQ'3 /( 1 - (f/fQ)2  + j2(f/fQ)} 

(1-7) 

Savage:  U3(f)  = fQ'3/{  1 - 3(f/fQ)2  + j[3(f/fQ)  - (f/fQ)3]}  . 

d-8) 

These  spectra  equal  f ^ at  0 frequency.  In  the  Sharpe  solution,  fQ  = a/(7r\T3a)  where  a is 
the  cavity  radius  and  a is  the  compressional  velocity  of  the  medium.*®  For  a constant  stress 
drop,  the  far-field  spectrum  scales  exactly  as  f **,  which  is  proportional  to  the  cavity  volume. 
This  geometric  similarity  also  is  assumed  to  hold  for  the  Brune  and  Savage  source  models. 

Figure  1-21  shows  the  amplitude  spectra  of  each  source  model  as  a function  of  scaled  fre- 
quency (f/fQ).  Sharpe  and  Brune  displacements  are  u;-square  models,  and  Savage’s  model  is 
a; -cube  at  high  frequency.  In  Fig.  1-22,  the  time-domain  waveforms  obtained  by  inverting 
Eqs.  (1-6),  (-7),  and  (-8)  are  displayed  as  a function  of  scaled  time  T = f t,  where  t is  real 
time.  These  displacement  waveforms  are  given  by 


Sharpe: 

u^(t)  = *JZvfo  exp  | j sin  (2^^  7rT) 

(1-9) 

Brune: 

u2(t>  = (2tt)2  fJ2T  exp  [— 2ttT] 

(1-10) 

Savage: 

u 3 (t ) = f^T2  exp  [ — 2 7r T]  . 

(I-U) 

The  area  of  each  waveform  in  normal  time  t equals  f the  spectral  level  at  0 Hz. 

Synthetic  records  were  generated  for  each  source  function  using  a set  of  corner  frequencies 

/ 1 

f = 0.0625,  0.125,  0.25,  . . . 16  Hz.  The  source  functions  were  convolved  with  the  Johnson- 
° 19 

Matheson  short-period  seismometer  response  used  at  the  Vela  arrays,  and  a causal  attenuation 

1 5 

operator  given  by  Carpenter.  Figure  1-23  shows  the  synthetic  waveforms  for  an  attenuation 
parameter  t*  =0.1,  normalized  to  the  same  gain,  t*  equals  travel  time  divided  by  average  Q. 

The  seismograms  for  the  Sharpe  and  Brune  source  functions  are  very  similar  with  a nar- 
row range  of  pulse  period^  from  ~0.3  to  0.8  sec.  On  the  other  hand,  the  seismograms  for  the 
Savage  source  show  a wider  range  of  periods  from  ~0.25  to  1.4  sec. 

Similar  synthetic  records  were  calculated  for  attenuation  parameters  t * = 0.25,  0.50,  and  1.0. 
Amplitudes  and  periods  were  measured  and  amplitudes  were  corrected  for  Vela  instrument 
response.  Corrected  values  of  log  (A/T)  are  displayed  in  Fig.  1-24  for  the  three  displacement 
models  as  a function  of  corner  frequency  f and  t*. 

Log  (A/T)  is  a measure  of  m^  except  for  an  overall  gain  factor.  In  Figs.  I-24(a)  and  (b),  a 

hatched  line  of  mH  vs  pulse  period  T for  U. S.  explosions  is  shown.  This  line  was  obtained  by 
20° 

P.  D.  Marshall  from  EDR  reports.  The  curves  of  log  (A/T)  of  the  Sharpe  model  were  arbi- 
trarily scaled  so  that  the  explosion  line  falls  in  the  region  between  the  curves  for  t * = 0.50  and 

i 

1.0,  a reasonable  assumption  based  on  direct  estimates  of  Q based  on  near-field  and  teleseismic 
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21  22 

measurements.  ' This  scaling  was  then  applied  to  the  theoretical  curves  in  Figs.  I-24(b) 
and  (c). 

It  is  apparent  that  the  range  of  periods  and  slope  of  the  explosion  line  matches  the  theoret- 
ical curves  for  the  Sharpe  and  Brune  models  reasonably  well.  However,  the  Rat  Island  earth- 
quake data  clearly  cannot  fit  either  co-square  model.  In  Fig.  I-24(c),  the  Savage  model  seems 
to  fit  the  earthquake  well  in  view  of  the  simplicity  of  the  model.  Both  the  wide  range  of  periods 
and  the  slope  of  m^  vs  T are  adequately  described  by  the  Savage  model.  Thus,  we  can  clearly 
reject  the  to -square  models  for  these  events  using  only  available  bulletin  data  and  without  spec- 
tral analysis.  C.W.  Frasier 

R.G.  North 


I.  CONTRIBUTION  OF  WWSSN  STATIONS  TO  THE  PDE  BULLETIN 

Three  separate  months  of  the  USGS  "Summary  of  Data  Entered  into  ERL  Hypocenter  System 
and  its  Usage"  reports  were  studied  to  evaluate  the  contributions  from  the  WWSSN  stations  used 
to  produce  the  Preliminary  Determination  of  Epicenters  (PDE)  Bulletin.  The  three  monthly  re- 
ports were  for  January  1972,  December  1972,  and  for  July  1973.  These  three  months  appear 
to  be  average  reporting  periods  to  the  USGS,  so  the  percentages  computed  of  the  total  events 
and  the  mean  day/event  of  reporting  times  could  be  used  as  the  approximate  percentage  and 
mean  day/event  time  for  any  given  period.  Each  station  of  the  WWSSN  network  was  evaluated 
as  a function  of  the  station’s  contribution  to  the  total  number  of  events,  its  contribution  to 
computations,  and  to  the  reporting  lag  time  of  the  station.  The  WWSSN  stations  used  for 
this  evaluation  are  listed  in  the  USGS  "Seismograph  Station  Abbreviations  and  Coordinates," 

March  1974,  deleting  stations  which  are  listed  as  closed. 

A total  of  115  WWSSN  stations  was  evaluated.  Thirty  of  these  stations  did  not  report  to 
the  USGS  in  any  of  the  above-mentioned  months.  The  percentage  of  the  stations'  time  picks 
associated  to  the  1548  events  reported  for  these  three  months  varied  from  43  percent  for  sta- 
tion KBL  to  0.3  percent  for  station  AKU,  with  the  network  station  average  association  being 
10.6  percent.  Although  many  of  the  WWSSN  stations  reported  amplitudes  and  periods  for  a 
large  percentage  of  their  reported  time  picks,  their  contributions  to  all  events  were  small  as 
the  number  of  associated  time  picks  was  small.  These  percentages  varied  from  13.3  percent 
for  station  COL  to  0.0  percent  for  28  stations.  The  average  WWSSN  station  contributed  ampli- 
tude and  period  measurements  for  3.7  percent  of  the  1548  events.  The  mean  time  in  days  for 
the  stations  of  the  WWSSN  network  to  report  events  varied  from  2.2  days  for  station  WES  to 
61.5  days  for  station  QUE.  The  mean  delay  for  the  entire  WWSSN  network  was  19.3  days.  Of 
the  8 5 stations  which  reported  during  these  three  months,  25  stations  reported  within  1 week, 

34  stations  within  2 weeks,  44  stations  within  3 weeks,  64  stations  within  4 weeks,  and  with  the 
remaining  21  stations  reporting  more  than  4 weeks  after  the  event. 

The  USGS  requires  that  at  least  5 station  arrival  times  be  associated  as  an  event  to  be  re- 
ported by  the  PDE  report.  A total  of  321  events  with  body  wave  magnitude  4.6  was  reported 
by  the  USGS  for  1972.  These  events  were  evaluated  to  determine  the  number  which  would  have 
met  the  5 associated  station  arrival  time  requirements  if  only  the  WWSSN  stations  were  reported. 
Two  hundred  fourteen  events  met  this  requirement  with  a loss  of  107  events  being  reported.  This 
does  not  mean  that  the  associations  of  W'WSSN  time  picks  for  these  214  events  would  meet  the 
travel-time  residual  restrictions  imposed  by  the  USGS  to  be  used  in  the  hypocenter  computation. 
One  hundred  fifty-three  events  of  the  321  events  would  have  met  the  required  number  of  stations 
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associated  as  well  as  the  restrictions  on  arrival  time  residuals.  Fifty-two  percent  of  the 
321  events  of  magnitude  4.6  would  not  have  been  reported  by  only  the  WWSSN  network  of 
stations.  Body  wave  magnitudes  for  107  of  these  153  events  would  have  been  calculated  from 
the  amplitudes  and  periods  reported  by  WWSSN  stations,  but  40  of  these  m^  values  would  have 
been  single -station  values. 

A further  investigation  to  determine  the  best  50  stations  using  existing  methods  of  report 
ing  to  the  USGS  is  being  conducted.  The  determining  factors  for  these  stations  will  be  the  num- 
ber of  times  the  arrival  times  of  these  stations  are  associated  and  used  to  compute  the  hypo- 
center,  the  number  of  amplitude  and  period  measurements  reported,  and  the  smallest  mean 
day/event  reporting  time.  R.  Needham 
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Fig.  1-1.  Mean  biases  measured  for  the  102  stations  used  in  this  study.  Circles  enlarge  as  biases 
become  more  negative,  crosses  enlarge  as  they  become  more  positive. 


Fig.  1-2.  Mean  biases  measured  for  stations  in  the  continental  U.S. 
Values  in  parentheses  are  from  Table  1-2,  all  others  from  Table  1-1. 
Dashed  lines  separate  regions  of  negative  and  positive  bias. 


Fig.  1-3.  Scatter,  as  defined  in  the  text,  of  station  m^  values  relative  to  average  m^, 
as  a function  of  average  m^,;  (a)  and  (b)  show  scatter  prior  to  and  after  application  of 
biases  as  station  corrections. 
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Fig.  1-4.  Seismic  source  regions  for  which  station  biases  have  been  calculated. 
The  number  of  events  in  each  reported  by  15  or  more  of  the  72  stations  in  Table  1-1 
is  shown  in  parentheses. 
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Fig.  1-5.  Annual  numbers  of  station  m^  reports  as  given  in  the  ISC  Bulletin 
for  (a)  all  stations,  (b)  72 -station  network,  and  (c)  13  stations  in  western  U.S. 


Fig.  1-6.  Ratios  of  annual  seismicity  as  a function  of  average  m^  for  1964-1969, 
1970-1973  to  that  for  1964-1973,  prior  to  and  after  application  of  biases  as 
station  corrections. 
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Fig.  1-7.  Plot  of  station  magnitude  bias 
(from  Table  1-1)  against  local  Pn  velocity 
(from  Marshall). 


Fig.  1-8.  Plot  of  station  magnitude  bias  (from  Table  1-1)  against  travel  time 
correction  (from  Herrin  and  Taggart^). 
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Fig.  1-9.  Magnitude  differences  for  stations  TFO  and  UBO,  using  UBO  as  a reference. 
Each  point  is  the  average  of  many  observations.  The  vertical  bars  show  95-percent 
confidence  bounds  on  the  mean  values.  The  solid  line  is  a theoretical  relation  derived 
from  parameters  explained  in  the  text. 


Fig.  1-10.  Magnitude  differences  for  stations  UBO  and  TFO,  using  TFO  as  a reference 
(compare  Fig.  1-9). 
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MAGNITUDE  (ALO)  - MAGNITUDE  (UBO) 


Fig.  I- 1 1 . Magnitude  differences  for  stations 
ALjQ  and  UBO,  using  UBO  as  a reference. 


Fig.  1-12.  Magnitude  differences  for  stations 
EUR  and  UBO,  using  UBO  as  a reference. 


28 


Fig.  1-13.  Corrected  values  of  log  (A/T) 
due  to  a Sharpe  source  model  as  a func- 
tion of  corner  frequency  and  attenuation 
parameter  t*.  The  dashed  line  is  the 
source  amplitude  at  1 Hz. 


Fig.  1-14.  Corrected  values  of  log  (A/T) 
due  to  a Brune  source  model  as  a func- 
tion of  corner  frequency  and  attenuation 
parameter  t*.  The  dashed  line  is  the 
source  amplitude  at  1 Hz. 


Fig.  1-15.  Corrected  values  of  log  (A/T) 
due  to  a Savage  source  model  as  a func- 
tion of  corner  frequency  and  attenuation 
parameter  t*.  The  dashed  line  is  the 
source  amplitude  at  1 Hz. 
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Fig.  1-16.  Bias  of  measured  log(A/T)  about  true  source  amplitude  at  1 Hz  for  Sharpe, 
Brune,  and  Savage  source  models.  Each  curve  shows  the  bias  for  a given  t*  as  a function 
of  corner  frequency.  The  m^  explosion  scale  was  obtained  from  Sharpe's  model  using 
empirical  relations  for  magnitude  m^,  yield  Y,  and  cavity  radius  a.  Numbers  in  paren- 
theses are  magnitudes  for  an  explosion  and  four  Aleutian  earthquakes  (see  text). 


30 


Fig.  1-17.  Map  centered  on  the  area  of  the  Rat  Island  sequence 
of  earthquakes  (1965).  Open  circles  are  WWSSN  stations  and  solid 
circles  are  Vela  arrays. 
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Fig.  1-18.  Amplitudes  and  periods  of  the  563  Rat  Island  earthquakes  recorded 
at  BMO:  (a)  shows  the  scatter  of  log  (A)  with  T;  (b)  displays  the  mean  and 

standard  deviation  of  log  (A)  at  each  T;  (c)  is  a histogram  of  reported  periods; 
and  (d)  is  a plot  of  mean  and  standard  deviation  of  log  (A/T)  as  a function  of  T. 


32 


log  (A/T) 


Fig.  1-19.  Mean  and  standard  deviation  of  log  (A/T)  as  a function  of  T 
for  Rat  Island  earthquakes  recorded  for  four  Vela  arrays. 
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Fig.  1-20.  Mean  and  standard  deviations  of  log  (A/T)  as  a function  of  T 
for  Rat  Island  earthquakes  recorded  at  four  WWSSN  stations. 
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Fig.  1-21.  Amplitude  spectra 
source  functions  plotted  as 
scaled  frequency  f/fQ,  where 
ner  frequency. 


Fig.  1-22.  Far-field  displacement  u(t)  for 
three  source  functions  plotted  as  a function 
of  scaled  time  tf  . 
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Fig.  1-23.  Synthetic  seismograms  obtained  by  convolving  three  source  functions 
with  the  Vela  short -period  seismometer  system.  Corner  frequencies  range 
from  fQ  = 0.0625  to  16  Hz.  Earth  attenuation  for  a t*  =0.10  is  included.  Am- 
plitudes of  all  seismograms  are  equalized  for  comparison. 
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T for  Rat  Island  earthquakes  recorded  simultaneously  at  all  five  Vela  stations. 


II.  SURFACE  WAVE  STUDIES 


A.  THE  EFFECTS  OF  LATERAL  HETEROGENEITY  IN  THE  OCEANIC  LITHOSPHERE 
UPON  SURFACE  WAVE  PROPAGATION 


Evidence  from  studies  of  Rayleigh  wave  dispersion  indicates  a systematic  increase  in  phase 

1 2 

velocity  as  the  age  of  the  sea  floor  increases  away  from  the  oceanic  ridge.  ' The  greatest 
change  takes  place  in  the  first  few  million  years,  due  to  the  rapid  cooling  and  solidification  of 

the  upper  part  of  the  lithosphere.  Anisotropy  of  propagation  such  that  the  Rayleigh  waves  travel 

1 

faster  in  the  direction  of  spreading  is  also  required  by  the  data  - this  feature  is  not  incorpo- 
rated here  since  it  considerably  complicates  the  differential  equations  to  be  solved  in  the  ray 
tracing  method  used. 

At  a period  of  3 3 sec,  the  values  obtained  for  Rayleigh  wave  phase  velocities  are  [for  age 
zones  0 to  5 million  years  (MY),  5 to  10  MY,  10  to  20  MY,  and  20  MY+]  3.794.  3.852,  3.920, 
and  3.966  km/sec,  respectively.  These  values  of  phase  velocity  may  be  assigned  to  the  mean 
ages  of  the  zones  (2.5,  7.5,  15,  and  25  MY),  where  the  last  is  a reasonable  value  since  very 
little  of  the  ocean  floor  studied^  is  older  than  38  MY.  These  data  points  are  shown  in  Fig.  II— 1. 

A very  good  fit  to  these  four  points  can  be  obtained  by  a second -order  polynomial  of  the  form 

V(A)  = aA2  + bA  + c . . . (II-l) 


where  V is  the  phase  velocity  at  age  A (MY).  Least-squares  values  of  (a,  b,  c)  are  (-0.239  X 

-3  -1 

10  , 0.1425  X 10  , 3.759).  The  value  of  phase  velocity  at  the  ridge  axis  (A  = 0)  is  given  by  a; 
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the  value  of  3.75  9 is  close  to  that  (3.769)  given  by  an  earth  model  for  the  mid -oceanic  ridge. 

The  ray  tracing  problem  is  considerably  simplified  in  this  context  if  we  consider  a spherical 
earth  composed  of  two  rigid  hemispherical  plates  such  that  their  axis  of  relative  rotation  coin- 
cides with  the  polar  axis  of  the  earth.  The  age  of  formation  A of  any  part  of  the  two  plates  at 
longitude  (p  is  then  given  by 


A = 


(II-2) 


where  <p^  is  the  longitude  of  the  ridge  axis  and  w is  the  relative  angular  rotation  vector  of  the 
two  plates.  The  age  A and,  thus,  phase  velocity  V(A)  are  then  conveniently  independent  of 
latitude  0.  The  relative  spreading  rate  at  the  ridge  axis  is  a maximum  at  the  equator  and  zero 
at  the  poles,  expressed  in  cm/yr;  it  is,  however,  constant  if  expressed  in  degrees  of  longitude 
<p/yr. 
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The  partial  differential  equations  for  surface  wave  ray  tracing  involve  0,  <p,  azimuth 
and  phase  velocity  V and  its  spatial  derivatives  8V/80  and  8V/8<p.  For  the  present  case,  taking 
(p^  = 0 and  a>  such  that  the  spreading  rate  is  1°  of  longitude/MY  (corresponding  to  11.1  cm/yr 
at  the  equator)  we  have,  from  Eq.  ( II - 1 ) 


V(cp)  = a (p 
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+ b (p  + c 
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( II— 3 ) 


making  the  ray  tracing  equations  particularly  simple. 
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The  results  of  ray  tracing  in  the  velocity  model  given  by  the  above  equations  are  shown  in 
Fig.  II-2.  The  rays  extend  to  a distance  of  25°,  corresponding  to  the  limit  of  validity  of  Eq.  (II-l). 
The  rays  leave  the  source  at  5°  intervals  in  azimuth  and  can  be  seen  to  be  slightly  bent  toward 
the  ridge  axis  (running  North-South).  The  velocity  derivative  bV/b<p  is  fairly  small  and  contin- 
uous; the  azimuth  derivation  of  the  ray  (azimuth  at  distance  of  25°  minus  starting  azimuth)  is  a 
maximum  (2.66°)  at  a starting  azimuth  of  8°. 

We  may  consider  the  error  introduced  in  determination  of  phase  velocity  as  a function  of 
age  (and  thus  longitude)  by  decomposition  of  the  path  into  great  circle  segments  as  opposed  to 
the  "true"  path  given  by  ray  tracing.  The  apparent  phase  velocity  at  a point  P is  given  by 


and 
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for  the  ray  path 
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for  the  great  circle  path 


/great  \ 
Vcircle/ 


(H-5) 


where  T(l)  is  the  phase  delay  introduced  by  a particular  segment  / of  the  path.  The  difference 
between  and  Yq^  points  P at  distances  of  2 5°  is  a maximum  at  the  azimuth  of  maximum 
azimuthal  deviation  of  8°;  it  does  not,  however,  exceed  0.0015  km/sec  and,  thus,  the  neglect 
of  horizontal  refraction  cannot  be  the  cause  of  the  apparent  anisotropy  observed.  Since  the  azi- 
muthal deviations  are  small,  the  effect  of  the  velocity  structure  upon  surface  wave  amplitudes 
is  also  small.  Figure  n-3  shows  the  variation  of  spectral  amplitudes  with  azimuth  due  to  non- 
geometrical  spreading.  The  amplitudes  can  be  seen  to  be  enhanced  in  the  direction  of  the  ridge 
axis.  This  implies  that  measurements  of  Q beneath  the  ridge  axis  from  surface  wave  ampli- 
tudes may  be  biased  upward  by  the  effect  of  nongeometrical  spreading. 


R.  G.  North 


B.  SYNTHETIC  RAYLEIGH  WAVE  SEISMOGRAMS  FOR  THE  FALLON  EARTHQUAKE 

The  event  pair  consisting  of  the  Fallon  earthquake  (20  July  1962)  and  the  Shoal  explosion 

(26  October  1963)  is  notable  because  it  is  one  of  the  few  occurrences  of  both  an  earthquake  and 

an  explosion  in  the  same  location  with  good  azimuthal  station  coverage.  Consequently,  it  has  a 

4-8 

long  history  of  inquiry  into  its  surface  wave  features.  This  present  study  of  the  Fallon  earth- 
quake is  an  attempt  to  derive  its  source  parameters  by  making  use  of  a seismic  moment  tensor 

q 

representation  of  the  Rayleigh  wave  point  source  equations.  Ultimately,  we  hope  to  obtain  these 
parameters  by  systematic  inversion.  The  results  presented  here  are  those  from  solutions  of 
the  forward  problem  (calculate  the  radiation  from  a given  moment  tensor)  and  were  found  by  an 
ad  hoc  trial  and  error  procedure. 

As  implemented,  the  theory  assumes  that  the  location  of  the  event  and  its  depth  are  known. 
Then,  using  a plane  layered  structure,  synthetic  Rayleigh  waves  are  calculated  for  the  case  of 
a step  function  of  applied  moment.  Other  source  time  functions  can  easily  be  accommodated  by 
the  appropriate  convolution  operation.  The  seismometer  instrument  response  is  also  included 
to  allow  direct  comparison  with  the  data.  The  synthetic  seismograms  presented  here  consist 
of  the  fundamental  and  first  higher  Rayleigh  wave  modes. 
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Both  Fallon  and  Shoal  occurred  in  the  Fairview-Dixie  Valley  region  of  Nevada  and  were  re- 
corded by  Long-Range  Seismic  Measurements  (LRSM)  stations  principally  situated  in  the  Basin 
and  Range  province.  Exceptions  are  the  California  stations.  Figure  II -4  shows  the  event  lo- 
cations and  the  11  LRSM  stations  which  were  used.  The  earth  structure  is  composed  of  a crust 

7 

appropriate  to  Fallon  (the  Fallon-Ruth  model  of  Toksoz  et  al.  ) underlaid  by  a typical  continen- 
tal mantle  with  both  a P-  and  S-wave  low  velocity  zone.  Phase  and  group  velocity  dispersion 
curves  as  a function  of  horizontal  wavenumber  (necessary  to  implement  the  theory)  are  shown 
in  Fig.  II-5.  The  data  plotted  to  the  same  amplitude  and  time  scales  are  shown  in  Fig.  II-6. 

Synthetic  seismograms  from  the  best  double -couple  fault  plane  solution  are  shown  in  Fig.  n-7. 
They  indicate  a double -couple  seismic  moment  near  1.1  X 10^  dyn-cm.  In  fitting  the  observed 
data,  it  was  noticed  that  the  solution  was  most  sensitive  to  the  slip  angle  A and  least  sensitive 
to  the  strike  angle  6.  Solutions  with  only  a few  degrees  different  slip  were  unacceptable,  but 
the  solution  persisted  over  15°  to  20°  of  strike  variation. 


TABLE  11-1 

FALLON  FAULT  PLANE  SOLUTIONS 

Strike 
(°E  of  N) 

Dip 

(deg  down) 

Slip 

(deg  up  HW  block) 

Depth 

(km) 

Present  study 

5 

75 

190 

10 

Lombert  et  ol . (1971) 

10 

82 

16 

15 

Toksoz  et  ol_.  (1965) 

355 

76 

230 

20 

A comparison  with  two  other  Fallon  fault  plane  solutions  is  presented  in  Table  II -1.  The 
present  study  solution,  in  addition  to  being  the  best  fit  to  the  data  in  the  vicinity  of  its  fault 

g 

angles,  comes  closest  to  matching  that  of  Lambert  et  al.  The  principal  differences  are  in  the 
direction  of  slip  (NNE  for  Lambert  et  aL,  SSW  for  the  present  study)  and  the  amount  of  dip  slip 
motion.  The  Fallon  data  given  here  favor  left  lateral  movement  with  relatively  little  dip  slip 

7 

motion.  The  solution  by  Toksoz  et  al.  is  based  on  only  three  stations  and  bears  little  resem- 
blance to  the  present  solution. 

There  is  some  tectonic  information  available  from  the  literature.  Van  Nostrand  and 
Helterbran^  show  surface  fault  breaks  from  a 1954  earthquake  in  the  Dixie  Valley  region  of 
Nevada.  The  Fallon  location,  computed  using  arrival  times  corrected  for  Shoal,  is  to  the  east 
Qf  one  of  the  1954  fault  breaks.  A crude  measurement  from  the  map  indicates  a strike  of  20° 
and  a dip  of  79°  which  tends  to  support  both  the  Lambert  et  al.  and  present  fault  plane  solutions. 

The  results  from  trying  to  add  compensated  linear  vector  dipole  (CLVD)  and  mono  pole  (MP) 
moment  sources  to  the  double-couple  (DC)  solution  are  shown  in  Fig.  II-8.  The  most  obvious 
effect  of  this  is  in  the  reduced  amplitude  of  the  WINV  seismogram.  The  addition  of  either  of 
these  source  mechanisms  to  the  DC  solution  tends  to  hurt  the  fit.  Thus,  the  actual  source 
mechanism  of  Fallon  must  have  appreciably  less  than  10  percent  CLVD  or  MP  components. 

These  results,  though  preliminary  in  that  they  represent  a subjective  solution  criterion, 
indicate  that  the  Fallon  earthquake  was  a left  lateral  predominantly  strike  slip  event  with  its 
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fault  plane  generally  aligned  with  known  tectonic  features  in  the  source  region.  The  double  - 

23 

couple  seismic  moment  was  approximately  1.1  X 10  dyn-cm  and  there  appeared  to  be  less 
than  10  percent  of  this  moment  in  either  the  CLVD  or  MP  mechanism.  Further  work  along 
this  line  should  consist  of  automating  the  inversion  process  partly  to  reduce  the  manual  labor 
involved  and  also  to  remove  the  subjective  judgment  factor. 

D.  W.  McCowan 


REFERENCES 


1.  D.  W.  Forsyth,  "The  Structural  Evolution  and  Anisotropy  of  the  Oceanic 
Upper  Mantle,"  Geophys.  J.  R.  Astr.  Soc.  43,  103-162  (1975). 

2.  A.  R.  Leeds,  "Lithospheric  Thickness  in  the  Western  Pacific,"  Phys. 
E^rth  Planet.  Inter.  11,  61-64  (1975). 

3.  B.  R.  Julian,  "Ray  Tracing  in  Arbitrarily  Heterogeneous  Media," 
Technical  Note  1970-45,  Lincoln  Laboratory,  M.  I.  T.  (31  December 
1970),  DDC  AD-720795. 

4.  J.  N.  Brune  and  P.  W.  Pomeroy,  "Surface  Wave  Radiation  Patterns 
for  Underground  Nuclear  Explosions  and  Small-Magnitude  Earth- 
quakes," J.  Geophys.  Res.  £8,  5005-5028  (1963). 

5.  "Fallon  Earthquake,"  Technical  Report  No.  63-122  (Rev.  1964),  Teledyne 
Geotech  Corporation,  Alexandria,  Virginia  (1964). 

6.  R.  Van  Nostrand  and  W.  Helterbran,  "A  Comparative  Study  of  the  Shoal 
Event,"  Seismic  Data  Laboratory  Report  109,  Teledyne  Geotech  Corpo- 
ration, Alexandria,  Virginia  (1964). 

7.  M.  N.  Toksoz,  D.  G.  Harkrider,  and  A.  Ben-Menahem,  "Determination 
of  Source  Parameters  by  Amplitude  Equalization  of  Seismic  Surface 
Waves,"  J.  Geophys.  Res.  70,  907-922  (1965). 

8.  D.  G.  Lambert,  E.  A.  Flinn,  and  C.  B.  Archambeau,  "A  Comparative 
Study  of  the  Elastic  Wave  Radiation  from  Earthquakes  and  Underground 
Explosions,"  Seismic  Data  Laboratory  Report  284,  Teledyne  Geotech 
Corporation,  Alexandria,  Virginia  (1971). 

9.  D.  W.  McCowan,  "Surface  Wave  Representation  of  Surface  WTave  Sources," 
Geophys.  J.  R.  Astr.  Soc.  44,  595-599  (1976). 


Fig.  II— 1.  Rayleigh  wave  phase  velocity  V 
at  33-sec  period  as  a function  of  age  of  for- 
mation A of  oceanic  lithosphere.  Points 
are  published  values  of  Forsyth  ^ for  2.5, 
7.5,  15,  and  25  MY;  line  is  polynomial  fit 

V(A)  = -0.239  X 10-3A2  + 0.1425  X 10"lA  + 
3.759. 
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RELATIVE  SPECTRAL  AMPLITUDE 


Fig.  n-2.  Result  of  ray  tracing  in  velocity 
model  given  by  velocity  function  of  previous 
figure,  for  spreading  rate  of  1°  longitude/ 
MY.  Polar  projection,  with  origin  centered 
on  ridge  axis:  great  circle  paths  are  straight 
lines  from  the  origin;  circles  denote  dis- 
tances of  5,  10,  15,  20,  and  25°. 


Fig.  II- 3.  Relative  spectral  amplitudes 
of  32- sec  Rayleigh  waves  as  a function 
of  azimuth  with  respect  to  ridge  axis, 
derived  from  Fig.  II -2  through  ampli- 
tude oc  (ray  density)!/^. 
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Fig.  II- 4.  Event  and  recording  station  locations  for  the  20  July  1962  Fallon  earthquake. 
Stations  are  temporary  LRSM  installations.  Data  set  consisted  of  long-period  vertical 
components. 
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Fig.  H-5.  Phase  and  group  velocity  dispersion  curves  for  72T  plane 
layered  earth  structure.  First  three  modes  shown.  Lines  of  equal 
period  are  hyperbolas. 
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Fig.  II-6.  Fallon  LPZ  data.  Station  codes  and  start  times  of  individual  traces 
are  given  at  the  left. 
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Fig.  H-8.  Synthetic  Rayleigh  wave  LPZ  seismograms  for  l6/l 8 DC,  l/l8  CLVD,  1/18MP 
source.  Plotted  to  same  scales  as  data.  Fault  angles  same  as  best  solution.  Total  moment 
is  1.1  X 10^3  dyn-cm.  Fundamental  plus  first  higher  modes  used  in  solution. 


III.  EVASION 


A.  DETECTING  MULTIPLE  EVENTS  IN  NARROWBAND  SEISMOGRAMS 

1 2 

It  has  been  suggested  * that  multiple  explosions,  suitably  timed,  can  be  used  to  simulate 
the  seismic  signature  of  an  earthquake.  Two  different  schemes  have  been  derived.  By  firing 
small  explosions  followed  by  several  large  ones  spaced  over  a few  seconds,  the  complexity  and 
Mg:m^  of  an  earthquake  can  be  mimicked  and  the  first  motion  obscured.  Alternatively,  the  ex- 
plosions may  be  spread  out  in  time  such  that  the  overall  trace  has  the  appearance  of  a first  ar- 
rival followed  by  a depth  phase.  By  suitable  timing,  the  second  arrival  can  even  be  made  to 
appear  inverted.  If  the  delay  time  indicates  a depth  deeper  than  is  technically  feasible  to  em- 
place an  explosive  device,  the  events  may  be  misidentified  as  a single  earthquake.  In  the  ab- 
sence of  close-in  stations  or  sufficiently  broadband  signal-to-noise  ratios  at  teleseismic  dis- 
3 

tances,  evasion  by  these  approaches  seems  possible  to  achieve  when  normal  discrimination 
techniques  are  employed. 

If  one  considers  that  in  both  instances  the  composite  seismograms  consist  of  approximately 
the  same  signal  delayed,  scaled,  and  summed  several  times,  the  problem  maybe  solved  by 
determining  the  arrival  time  of  each  signal,  and  so  showing  the  trace  to  be  the  result  of  multi- 
ple events  rather  than  a single  event.  What  we  wish  to  do,  then,  is  to  deconvolve  observed 
seismograms.  If  the  source  function  is  known,  the  problem  reduces  to  the  classical  Wiener 
filtering  problem.  Unfortunately,  signal  shapes  from  various  explosions  are  too  different  for 
that  technique  to  be  definitive.  When  neither  the  source  nor  the  echo  string  is  known,  the  prob- 
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lem  is  referred  to  as  blind  deconvolution.  The  solution  to  such  a problem  can  be  obtained 
through  homomorphic  deconvolution.5  When  the  source  function  is  long  compared  to  the  size  of 
the  echo  train,  many  segments  of  log  spectral  estimates  may  be  averaged  to  yield  the  log  spec- 
trum of  the  echo  train  alone.  For  the  previously  mentioned  evasion  scenarios,  the  echo  train 
and  the  source  are  comparable  in  size.  In  addition,  the  recording  response  of  short-period 
seismometers  is  narrow  and  the  log  spectrum  is  uncontaminated  by  noise  only  over  that  small 
interval.  In  such  a case  maximum  entropy  cepstral  analysis  ’ is  an  appropriate  method  to  ef- 
fect the  deconvolution.  The  following  experiments  were  conducted  to  illustrate  the  usefulness 
of  using  maximum  entropy  cepstral  analysis  on  narrowband  short-period  teleseismic  record- 
ings to  detect  multiple  arrivals  that  appear  as  a single  earthquake  recording. 

Figure  III-l(a)  shows  the  recording  of  a single  Nevada  Test  Site  (NTS)  explosion  observed 

at  NORSAR.  Figure  III-l(b)  shows  the  complexity  simulation  and  Fig.  III-l(c)  the  depth  phase 

2 

simulation,  according  to  the  delays  given  by  Landers.  If  the  seismogram  in  Fig.  Ill-  1(a)  is 
denoted  by  S(t),  then  the  signature  in  Fig.  Ill- 1 (b ),  C(t)  is 

C(t)  = S(t)  © [26(t)  + 56 (t  - 1)  + 9 6 (t  + 2) 

+ 206 (t  - 4)  + 206 (t  - 6)  + 406(t  - 7)  + 306(t  - 8)  + 256 (t  - 9)] 

and  in  Fig.  Ill- 1 (c),  D(t)  is 

D(t)  = S(t)  © [86(t)  + 3 6 (t  - 7)  + 7 6 (t  - 7.7)  + 46(t  - 8.1) 

+ 56(t  - 10.5)] 

where  © signifies  convolution. 

For  C(t),  the  record  is  sufficiently  complex  so  that  such  a seismogram  could  not  have  been 
written  if  the  source  were  a single  explosion.  Additionally,  the  body  wave  magnitude,  as 
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conventionally  measured  in  the  first  cycle,  would  be  approximately  0.6  smaller  than  that  due  to 
the  largest  single  explosion.  The  surface  wave  magnitude  would  be  increased  by  approximately 
0.4  over  that  for  the  largest  single  explosion,  since  the  energy  from  all  the  explosions  would 
simply  sum  in  the  20-sec  band.  The  combination  of  the  decrease  in  and  increase  in  Mg 
would  result  in  their  ratio  falling  in  the  earthquake  population.  Taking  the  log  spectrum  of  C(t), 
correcting  for  instrument  response,  and  using  maximum  entropy  spectral  analysis  (in  the  good 
signal -to- noise  band,  0.2  to  4.5  sec),  the  cepstrum  of  C(t)  shown  in  Fig.  Ill- 2 is  obtained.  The 
1-  and  2-sec  delays  clearly  indicate  the  multiple  arrivals.  The  depth  phase  in  S(t)  is  near  1 sec 
and,  so,  not  resolvable. 

For  the  depth  phase  simulation,  D(t)  appears  to  be  composed  of  the  P,  pP  and  sP  phases 
for  a single  event  whose  depth  is  near  3 0 km,  thus  indicating  that  the  event  is  an  earthquake. 
Using  that  band  of  D(t)  near  the  apparent  depth  phase  pP,  the  cepstrum  shown  in  Fig.  Ill— 3 is 
obtained.  The  cepstrum  indicates  a triple  event  consisting  of  a smaller  arrival  preceding  the 
major  arrival  by  0.7  sec  and  a larger  arrival  following  it  at  0.5  sec.  These  delays  are  in  good 
agreement  with  the  actual  delays  of  0.7  and  0.4  sec,  respectively. 

We  conclude,  then,  that  for  explosions  at  the  lower  magnitudes  where  the  usefulness  of 
broadband  recording  is  negated  by  poor  signal -to- noise  ratios,  the  occurrence  of  multiple  events 
appearing  as  single  earthquakes  on  teleseismic  recordings  can  be  deciphered.  The  techniques 
used  in  these  experiments  and  the  advent  of  the  Seismic  Research  Observatory  (SRO)  network 
from  which  the  necessary  high-quality  digital  data  can  be  obtained,  make  evasion  by  multiple 
events  a difficult  goal  to  achieve.  T y ^ anders 


B.  COMPARISON  OF  DIFFERENT  PREDICTION  ERROR  FILTERS 
ON  A SYNTHETIC  COMPOSITE  EVENT 


In  an  attempt  to  determine  the  effectiveness  of  prediction  error  filtering  in  unraveling  com- 
posite events,  we  applied  different  types  of  filters  to  a synthetic  event  made  up  of  multiple  re- 
cordings of  an  Eastern  Kazakh  presumed  explosion.  The  presumed  explosion,  recorded  by  a 
LASA  D2  subarray  SPZ  instrument,  was  on  24  July  1970  and  is  shown  as  the  top  trace  in 
Fig.  III-4. 

The  first  part  of  the  experiment  consisted  of  designing  a 30-point  (li  sec)  unit  span  predic- 
tion error  filter  using  3 0 sec  of  this  waveform.  The  algorithm  was  the  usual  least-squares 
Wiener-Hopf  method  which  has  been  successfully  used  in  the  petroleum  industry  for  many  years. 
The  output  from  this  filter  is  shown  as  the  bottom  trace  in  Fig.  III-4.  The  filter  condenses  the 
original  waveform  which  is  approximately  4 sec  long  to  a compact  waveform  consisting  of  a 
small  positive  spike  followed  by  a larger  negative  spike  0.33  sec  delayed  followed  by  a broad 
positive  excursion  0.76  sec  delayed.  The  total  length  of  the  condensed  waveform  is  on  the  order 
of  1 sec. 

The  synthetic  composite  event  was  constructed  using  the  amplitudes  and  delays  formulated 
2 

by  Landers.  These  were  selected  to  simulate  the  pP  phase  in  deep  events  and  are  given  below: 


Phase  Relative  Amplitude  Delay  (sec) 


1 8 0.0 

2 3 4.0 

3 7 4.7 

4 4 5.1 

5 5 7.5 
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The  synthetic  composite  event  (first  trace)  and  the  results  of  using  30-point  unit  span  prediction 
error  filters  designed  on  30  sec  of  the  original  waveform  (second  trace)  and  30  sec  of  the  com- 
posite waveform  (third  trace)  are  shown  in  Fig.  ITI-5. 

As  can  be  seen,  there  is  very  little  difference  in  the  results  between  designing  the  filters 
on  the  original  vs  the  composite  waveform.  In  general,  however,  the  results  will  not  be  this 
similar  when  an  actual  composite  event  is  used  because  the  noise  will  not  have  been  delayed 
and  summed  the  same  way  the  composite  was  constructed.  By  keying  on  the  large  negative  spike, 
four  of  the  five  phases  used  to  construct  the  composite  can  be  readily  identified  and  their  time 
delays  measured  to  within  5 percent.  To  eliminate  the  subjective  judgment  element  as  much  as 
possible,  these  results  were  submitted  to  a trained  analyst  (L.  Lande)  for  examination.  He,  too, 
was  able  to  identify  and  measure  four  of  the  five  phases  in  the  filtered  output.  The  fifth  phase 
(the  fourth  in  the  tabulation  above)  can  be  seen  in  the  output  only  when  its  time  delay  is  known. 
Since  this  phase  occurs  with  the  smallest  delay  between  itself  and  its  predecessor  (0.4  sec),  we 
conclude  that  such  a small  delay  is  beyond  the  resolution  of  our  method. 

The  composite  event  was  then  processed  using  a 30-point  prediction  error  filter  designed 

7 

with  the  adaptive  algorithm  described  by  Griffiths  et  al.  The  results  for  three  different  learn- 
ing rates  are  shown  in  Fig.  Ill- 6.  As  shown  by  Griffiths,  small  values  of  a (<  0.5)  produce  slow 
learning  rates  but  good  convergence  in  the  limit  and  large  values  of  a (<  1.0)  produce  rapid  but 
noisy  convergence  to  the  optimum  least-squares  filter.  Again,  four  out  of  the  five  phases  in  the 
composite  can  be  identified  and  their  time  delays  accurately  measured.  The  best  results,  in 
terms  of  balance  between  noise  and  resolution,  appear  to  be  with  a = 0.5. 

To  summarize,  we  have  shown  from  these  results  in  prediction  error  filtering  a synthetic 
composite  event  that  either  the  least  squares  or  adaptive  algorithms  is  able  to  provide  the  reso- 
lution necessary  to  identify  four  out  of  the  five  individual  phases  in  our  composite.  Neither 
method  can  resolve  the  last  phase  occurring  at  a relative  delay  of  0.4  sec  to  its  predecessor. 
Furthermore,  the  results  between  designing  a least-squares  filter  on  the  original  vs  the  com- 
posite waveforms  are  quite  similar.  D McCowan 
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Fig.  Ill— 1 . Basic  source  function, 
complexity  simulation,  and  depth 
phase  simulation,  respectively, 
discussed  in  the  text:  (a)  S(t), 

(b)  C(t),  and  (c)  D(t). 
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Fig.  Ill  — Z . Maximum  entropy 
cepstrum  of  C(t). 


Fig.  Ill  — 3.  Maximum  entropy  cepstrum 
of  the  apparent  depth  phase  pP  portion 
of  D(t). 
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Fig.  Ill -4.  Presumed  explosion  from  Eastern  Kazakh,  24  July  1970,  recorded 
at  the  D2  LASA  subarray  prediction  error  filtered  with  a 30 -point  unit  span 
predictor  designed  on  30  sec  of  the  waveform.  Five  percent  white  noise  has 
been  added  to  the  filter  design  equations. 


PREDICTION  ERROR  FILTER 
DESIGNED  ON  ORIGINAL 
WAVEFORM 


PREDICTION  ERROR  FILTER 
DESIGNED  ON  COMPOSITE 


Fig.  Ill  - 5.  Composite  event  prediction  error  filtered  with  30 -point  unit  span  predictors 
designed  on  both  the  original  waveform  and  the  composite  event  waveform.  Five  percent 
white  noise  has  been  added  to  the  filter  design  equations. 
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Fig.  Ill -6.  Composite  event  prediction  error  filtered  with  30 -point  unit  span  adaptive 
predictors.  Constant  a controls  the  learning  rate  of  the  predictor. 
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IV.  EARTH  HETEROGENEITY 


A.  ANALYSIS  OF  SCATTERING  FROM  NOVAYA  ZEMLYA  EXPLOSIONS 

1 2 

In  two  previous  SATS,  ' we  studied  the  codas  of  three  component  short-period  data  from 
six  presumed  explosions  at  Novaya  Zemlya  recorded  at  LASA.  Using  a polarization  filter,  it 
was  shown  that  P codas  are  extremely  sensitive  to  small  changes  in  source  location,  implying 
that  observed  body  phase  components  are  not  primarily  generated  by  structure  under  LASA. 

We  have  continued  this  study  by  investigating  the  arriving  azimuth  and  emergence  angle  of 
the  strongly  polarized  phases  in  the  codas,  in  an  attempt  to  isolate  possible  scattering  regions 
in  the  earth.  First,  we  pick  the  strongest  arrivals  in  the  coda  by  calculating  the  amplitude  of 
the  particle  motion.  An  envelope  is  constructed  by  convolving  the  instantaneous  amplitude  of 
the  polarized  waves  with  a triangular  window.  The  envelope  is  then  fitted  with  a least-squares 
exponential  decay  to  approximate  the  average  coda  level  with  time.  An  example  of  this  is  shown 
in  Fig.  IV-1  using  event  1,  which  occurred  on  14  October  1969  with  an  m^  of  6.1.  All  envelope 
peaks  which  are  above  a fraction  of  the  exponential  decay  curve  plus  background  noise  are  se- 
lected. The  apparent  dip  and  azimuth  of  the  selected  peaks  are  then  calculated,  and  these  angles 
are  plotted  on  an  equal  area  projection.  Figures  IV-2  to  -5  are  some  examples  of  these  pro- 
jections for  event  1,  each  figure  representing  a minute  of  data  after  the  P arrival.  Figure  IV-6 
shows  an  equal  area  projection  for  selected  amplitude  peaks  in  the  fourth  minute. 

If  the  coda  waves  are  generated  near  source  region,  the  emergence  angles  should  be  con- 
centrated around  the  initial  P-wave  direction.  Some  projections,  however,  show  that  the  waves 
arrive  from  all  directions,  and  many  phases  have  very  shallow  apparent  emergence  angles  (see 
Fig.  IV -5).  For  a P-wave  incident  on  a free  surface,  due  to  the  reflected  SV  wave,  the  apparent 
emergence  angle  will  be  about  19°  even  at  near  zero  true  emergence  angle  for  a medium  obey- 
ing Poisson's  relation.  Thus,  the  peaks  showing  up  near  the  rim  of  the  equal  area  projection 
are  most  likely  to  be  shear  waves. 

The  existence  of  these  shear  waves  shows  that  a considerable  part  of  the  coda  energy  is  not 
direct  P or  P-to-P  scattering  from  regions  near  the  source.  The  arrival  time  of  these  shear 
waves  in  coda  is  not  consistent  with  any  known  shear  phases  in  travel-time  tables.  It  is  unlikely 
that  the  shear  phases  in  Figs.  IV-5  and  -6  are  P-to-S  conversion  under  LASA,  because  there  are 
not  so  many  P phases  preceding  them  in  Fig.  IV-4.  It  is  thus  likely  they  are  scattered  shear 
waves  of  P-to-S  type  from  somewhere  along  the  path. 

Statistically,  shear  energy  is  important  for  scattering  of  P-waves  by  small  inhomogeneities 
only  when  ka  is  not  much  larger  than  one,  where  k is  the  wavenumber  and  a is  the  correlation 
distance  of  inhomogeneities.  ' If  we  assume  the  P-wave  velocity  is  in  the  order  of  6 km/sec, 
the  existence  of  scattered  S-wave  about  1 Hz  would  mean  that  somewhere  along  the  P-wave  hit 
a region  of  inhomogeneity  with  correlation  distance  not  much  larger  than  1 km. 

It  is  thus  inferred  that  the  impulsive  looking  body  phases  in  codas  are  scattered  from  in- 
homogeneities in  a rather  broad  region  not  limited  to  near  source  or  near  receiver  regions. 

The  short  duration  and  high  linearity  of  these  phases  indicate  that  the  area  of  these  inhomogene- 
ities is  not  too  large. 

Lacking  an  array  with  three  component  short-period  seismometers,  at  the  present  time,  we 
cannot  locate  the  exact  scattering  region  for  these  phases.  However,  we  can  try  to  match  strong 
coda  arrivals  with  broad  scattering  regions  in  the  earth  using  travel-time  tables.  To  illustrate 
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this  we  chose  two  phases  which  appear  on  the  records  of  all  four  northern  events.  In  Fig.  IV -7, 
a possible  P-to-P  scattered  phase  arrives  at  1 min.  17  sec  after  P,  and  in  Fig.  IV-8,  a possible 
P-to-S  scattered  phase  occurs  at  2 min.  55  sec  after  P.  Figures  IV-9  and  -10  show  the  epi- 
central  lines  of  scattering  regions  which  produce  P-to-P  and  P-to-S  arrivals  at  the  times  given 
above.  Each  line  is  a projection  onto  the  earth's  surface  of  the  scattering  locus  at  the  indicated 
depth.  It  can  be  seen  that  the  compatible  S-scattering  region  is  limited  within  20°  from  LASA, 
while  the  compatible  P-scattering  region  extends  to  larger  distances. 

M.  Yang 

C.  W.  F rasier 


B.  MISLOC ATION  PATTERNS  FOR  SEISMIC  NETWORKS 

In  previous  communications,  attention  has  been  drawn  to  the  array  diagram  at  LASA^  and 
NORSAr/  Lateral  heterogeneity  in  the  deep  mantle  and  near  the  source  region  or,  alternatively, 
possibly  small-scale  inhomogeneities  in  the  upper  mantle  under  the  arrays/  have  been  proposed 
to  explain  certain  anomalous  features  in  the  pattern  of  slowness  vectors.  Since  the  effect  of 
small-scale  upper-mantle  inhomogeneities  diminishes  with  increasing  array  aperture,  it  makes 
sense  to  compare  conventional  array  data  with  data  from  an  extended  network  of  stations  in  the 
same  region.  Here  we  compare  mislocations  for  NORSAR  with  those  for  the  Scandinavian  net- 
work, and  mislocations  for  LASA  with  those  for  a network  in  Montana.  Stations  reporting  to 
ISC  have  been  used;  the  Scandinavian  network  also  includes  Finland;  the  Montana  network  has 
been  restricted  to  latitudes  38°  to  54°N,  longitudes  98°  to  11  5° W.  Reference  locations  are  from 
NOAA  and  arrival  time  data  from  the  ISC  bulletins. 

Previously,  array  mislocations  have  been  obtained  from  slowness  vector  anomalies/  With 
increasing  aperture  of  the  array  or  network,  sphericity  of  the  earth's  surface  and  curvature  of 
the  wave  front  become  important.  These  aspects  may  be  taken  into  account  by  computing  slow- 
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ness  vectors  from  a parameterized  travel  time  response.  Alternatively,  mislocations  may  be 
obtained  from  a first-order  perturbation  of  the  initial  event  latitude  and  longitude.  Provided  a 
suitably  damped  iteration  scheme  is  used,  the  latter  method  more  often  converges  to  a stable 
solution.  Stability  criteria  were:  RMS  time  residue  <1.4  sec,  number  of  stations  >5.  Typi- 
cally, number  of  stations  were  10  in  Scandinavia  and  7 in  Montana.  This  method  has  been  used 
and  applied  to  three  source  regions  (China-Russia  border,  Central  America,  Bonin  Arc  and 
surroundings),  where  pronounced  and  rapidly  varying  slowness  anomalies  at  LASA  and  NORSAR 
have  been  reported/*6  Figures  IV-11  to  -13  summarize  the  resulting  mislocations,  which  may 
be  compared  with  similar  previously  published  figures  from  LASA  and  NORSAR  data.  Some 
points  in  this  comparison  are  noted  here: 

(1)  China-Russia  Border  Region:-  From  the  Hindu  Kush  to  the  explosion  site 
in  Eastern  Kazakh,  about  1700  km  to  the  North-East  and  corresponding 
with  only  a few  degrees  in  azimuth  from  NORSAR,  there  is  a complete  re- 
versal of  the  large  NORSAR  mislocations/  Data  from  the  Scandinavian 
network  (Fig.  IV-11)  show  no  significant  anomalies  throughout  the  entire 
region  considered.  A preliminary  evaluation  of  this  result  places  the 
source  of  the  NORSAR  anomaly  in  the  near  array  structure. 
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(2)  Central  America:-  NORSAR  mislocations  vary  significantly  with  geo- 
graphic region  throughout  Central  America.  In  particular,  a pro- 
nounced splitting  in  location  occurs  for  events  in  the  Oaxaca  district  of 
Mexico  (near  15°N,  90 °W).^  The  overall  trend  of  mislocations  in  the 
same  region  obtained  with  the  Scandinavian  network  (Fig.  IV -12)  appears 
to  be  radically  different,  but  some  important  details  are  remarkably 
similar  (discarding  6 unusual  large  arrows,  which  correspond  to 
dT/dA  near  or  below  the  core  shadow  boundary  value).  For  example, 
there  is  clear  evidence  for  lateral  variation  in  dT/dA  anomalies  and, 

in  particular,  the  splitting  of  NORSAR  mislocations  near  (15°N,  90 °W) 
also  seems  to  be  present  in  Fig.  IV-12.  This  result  supports  the  sug- 
gestion of  a related  anomaly  in  the  deep  mantle  on  the  source  side  of 
the  ray  path. 

(3)  Bonin  Arc  Region:-  LASA  mislocations  of  events  in  the  Bonin  Arc  region 
contrast  sharply  with  mislocations  of  neighboring  events  in  the  regions 

of  the  Ryukyu  and  Mariana  Islands.^  Mislocations  with  the  Montana  net- 
work (Fig.  IV -13)  do  not  reflect  a proposed  source  related  Bonin  anomaly. 
There  are  other,  less  conspicuous,  features  which  need  further  analysis. 
For  example,  locations  of  Mariana  Islands  events  with  the  Scandinavian 
network  (not  shown  here)  exhibit  a splitting  phenomenon  similar  to,  though 
less  striking  than,  the  Central  American  anomaly. 


In  summary,  most  of  the  large  array  anomalies  seem  to  be  related  to  structure  near  the 
receiver,  but  some  have  their  cause  in  the  deep  mantle  or  near  the  source  region.  The  latter 
anomalies  were  surprisingly  well  resolved  by  the  localized  networks,  in  view  of  the  scatter  in 
bulletin  data,  variability  in  effective  number  of  stations,  and  the  lack  of  suitable  station  correc- 
tions. It  is  to  be  expected  that  more  interesting  features  emerge  from  other  source -receiver 


combinations. 


D.J.  Doornbost 
J.  M.  Vermeulem 


C.  LARGE-SCALE  HETEROGENEITIES  IN  THE  LOWER  MANTLE; 

CORRELATION  WITH  THE  GRAVITY  FIELD 

In  our  previous  report,1  we  described  preliminary  results  of  a three-dimensional  inversion 
of  over  700,000  P-wave  travel-time  residuals  taken  from  the  ISC  Bulletin  for  the  years  1964- 
1970. 

During  the  last  year,  we  have  reanalyzed  the  data  using  different  discretization  patterns. 
The  shell  boundaries  were  changed  with  respect  to  those  listed  in  Table  IV -1  of  Ref.  1,  by  sub- 
division of  Region  II  into  two  shells:  Region  IIA  extending  from  670-  to  11 00 -km  depth  and  Re- 
gion IIB  from  1100  to  1500  km.  The  upper  boundary  of  Region  III  was  lowered  from  1400-  to 
1500 -km  depth.  This  subdivision  of  Region  11  has  significantly  reduced  the  compensation  in 
velocity  anomalies,  obvious  in  Figs.  IV-l(a)  and  (b)  of  Ref.  1.  It  appears  that  there  is  an  im- 
portant change  in  the  pattern  of  velocity  anomalies  at  a depth  of  approximately  1000  km.  With 


t Vening  Meinesz  Laboratory  for  Geophysics  and  Geochemistry,  University  of  Utrecht, 
The  Netherlands. 
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the  horizontal  boundaries  located  as  in  Ref.  1,  we  have  now  150  unknown  parameters  6v;  wc 

shall  identify  the  result  of  this  inversion  as  Case  I. 

In  order  to  test  the  stability  of  the  solution  with  respect  to  a change  in  location  of  the  bound* 

aries,  we  rotate  the  meridional  boundaries  by  30°;  the  result  of  this  inversion  is  identified  as 

o 

Case  II  (we  do  not  present  the  full  numerical  details  here;  the  complete  report  is  now  in  the 
process  of  publication). 

From  the  beginning,  we  have  intended  to  represent  lateral  variations  of  velocity  in  terms 
of  the  spherical  harmonic  expansion;  the  discretization  approach  has  been  adopted  only  for 
reasons  of  computational  efficiency. 

In  the  process  of  inversion,  we  have  determined  velocity  perturbations  <5v.^,  with  standard 
errors  a.  where  i is  the  shell  index,  and  j and  k are  the  indices  of  latitudinal  and  longitu- 
dinal zones,  respectively.  If  all  the  <5v  should  have  equal  weight,  the  simplest  way  of  deter- 
mining coefficients  of  the  spherical  harmonic  expansion  is  to  evaluate  integrals: 


A.  t 
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We  call  the  set  of  coefficients  obtained  in  this  way  the  T solution.  It  is  also  possible  to 
obtain  coefficients  A and  B by  the  least-squares  solution  of  an  overdetermined  system  of  ob- 
servational equations: 
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weighted  by  a.^.  The  bars  over  the  P^  and  cos  m<p  or  sin  nup  designate  the  average  value  for 
the  j**1  latitudinal  or  kth  longitudinal  zone.  It  would  appear  that  the  weighted  least-squares  so- 
lution (W)  should  be  preferred  over  the  T solution,  but  the  values  of  a.^  must  be  treated  with 
some  caution,  as  the  sampling  of  individual  blocks  may  be  uneven.  Thus,  differences  between 
the  sets  of  A.^m  and  obtained  using  Eqs.  (IV-1)  and  (IV-2)  represent  a measure  of  uncer- 

tainty in  the  actual  values  of  the  spherical  harmonic  coefficients. 

We  compute  cross-correlation  coefficients: 
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(1)  between  the  T and  W solutions  for  each  case  separately  to  establish  the  level  of  consistency 
of  the  velocity  perturbations  depending  on  the  method  of  derivation  of  coefficients,  and  (2)  be- 
tween Case  I and  II  coefficients,  separately  for  the  T and  W solutions,  to  investigate  the  sta- 
bility of  the  solution  with  respect  to  a change  in  the  discretization  of  the  model.  The  coeffi- 
0 3 3 

cients  A^  are  excluded;  the  coefficients  B^  and  A^  are  excluded  when  cross-correlation  between 
Case  I and  II  solutions  is  evaluated,  as  one  is  indeterminate  for  each  solution,  respectively. 

In  Table  IV-1,  we  list  for  these  four  cases  the  diagonal  elements  of  the  cross-correlation 
matrix.  We  note  that  in  most  cases  the  cross -correlation  coefficients  between  the  T and  W 
solutions  are  on  the  order  of  0.9;  this  means  that  these  two  different  methods  do  not  lead  to 
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TABLE  IV- 1 

CROSS-CORRELATION  COEFFICIENTS  BETWEEN 
T AND  W SOLUTIONS 

1 T/l  Wt 

II  T/ll  W 

1 T/ll  T 

1 W/ll  w 

Region  1 

0.960 

0.793 

0.085 

0.276 

Region  IIA 

0.913 

0.813 

0.064 

0.009 

Region  IIB 

0.854 

0.502 

0.315 

0.541 

Region  III 

0.881 

0.833 

0.895 

0.930 

Region  IV 

0.875 

0.855 

0.895 

0.828 

t Cross- correlation  between  T solution,  Cose  1,  and  W solution. 
Case  1.  Other  columns  are  labeled  similarly. 

grossly  different  representations  of  velocity  anomalies  in  terms  of  spherical  harmonics.  We 
also  note  that  we  obtain  coefficients  of  cross -correlation  between  the  Case  I and  Case  II  solu- 
tions that  are  nearly  zero  for  Regions  I and  IIA,  assume  intermediate  values  for  Region  IIB, 
and  are  high  (0.9)  for  Regions  III  and  IV. 

A plausible,  but  by  no  means  unique,  interpretation  of  this  pattern  of  correlation  coeffi- 
cients is  that  lateral  heterogeneities  in  the  mantle  above  1100  km  depth  are  dominated  by  short 
wavelength  features,  with  dimensions  considerably  below  those  of  our  grid  pattern;  the  inco- 
herent results  represent  the  effect  of  a wavenumber  bias  due  to  uneven  sampling  by  sources 
and  receivers.  The  mantle  below  1500  km  on  the  other  hand,  could  be  dominated  by  heterogene- 
ities of  dimensions  comparable  with  our  block  size,  or  greater,  as  the  results  do  not  appear  to 
be  sensitive  to  the  position  of  boundaries  between  blocks.  There  may  be  a transition  in  the  scale 
of  heterogeneities  in  the  depth  range  between  1100  and  1500  km. 

In  order  to  derive  the  final  T and  W solutions,  we  combine  velocity  perturbations  deter- 
mined for  Case  I and  Case  II.  For  the  integration  method  (T),  this  corresponds  simply  to  aver- 

3 3 

aging  the  coefficients  obtained  for  both  grid  patterns,  except  for  the  terms  and  . The 
weighted  least- squares  solution  (W)  is  obtained  by  solving  the  system  of  60  observational  equa- 
tions for  each  layer.  The  spherical  harmonic  coefficients  obtained  in  this  way  are  listed  in 
Table  IV -2. 

The  decay  with  angular  order  number  of  the  magnitude  of  the  spherical  harmonic  coefficients 
of  the  observed  gravitational  potential  suggests  that  the  source  of  the  low-order  anomalies  is 
located  in  the  deep  mantle;  such  an  interpretation  is,  of  course,  not  the  only  one  possible. 

We  test  this  hypothesis  by  assuming  that  the  density  perturbations  should  be  proportional  to 
velocity  perturbations.  We  supplement  Eq.  (IV-2)  by  equations  of  constraint: 

£HilAilm  = “C<m  ; 


= QS, 

it  lim  im 


(IV-4) 
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TABLE  IV-2 

SPHERICAL  HARMONIC  EXPANSION  COEFFICIENTS  OF  FINAL  T AND  W SOLUTIONS 


Region  1 

Region  1 1 A 

Region  MB 

Region  III 

Region  IV 

T 

W 

T 

w 

T 

W 

Gt 

T 

W 

Gt 

T 

W 

Gt 

O o 
< 

1.9 

1.6 

-3.3 

-1.5 

22.0 

21.0 

23.9 

-9.8 

-10.3 

-8.5 

24.7 

22.2 

24.1 

1.2 

1.1 

-1.5 

-2.6 

1.6 

3.5 

-0.3 

-0.4 

0.5 

-1.7 

4.2 

6.3 

4.7 

*; 

-1.3 

-1.8 

2.3 

1.2 

0.4 

0.0 

-1.1 

-0.8 

-0.6 

-0.8 

-12.8 

-12.3 

-11.4 

-1.5 

-0.8 

2.1 

-0.1 

0.4 

1.2 

0.2 

0.3 

0.2 

-0.1 

-1.9 

-2.3 

-1.8 

O CN 
< 

-0.4 

0.1 

0.9 

-0.4 

2.4 

0.6 

2.9 

5.2 

4.5 

5.8 

0.6 

-0.1 

1.7 

-0.7 

-0.9 

-2.0 

0.3 

1.0 

4.9 

5 .9 

-0.4 

-0.5 

0.1 

4.9 

1.4 

2.3 

1.6 

3.8 

-0.4 

-1.3 

-2.3 

-2.3 

0.0 

0.9 

1.2 

3.2 

1.8 

1.8 

3.1 

A2 

0.6 

-1.6 

2.3 

7.6 

-8.6 

-9.7 

-10.6 

2.0 

4.8 

5.0 

-0.8 

-2.5 

-0.6 

-0.4 

-0.8 

1.3 

2.2 

-1.7 

-1.6 

3.0 

-2.9 

-4.0 

-1.3 

-1.7 

-0.1 

2.7 

*3° 

3.1 

4.6 

0.3 

0.6 

0.2 

1.0 

2.3 

-4.9 

-3.5 

-2.6 

-6.1 

-8.5 

-8.0 

A3 

1.8 

3.1 

0.7 

-1.8 

3.2 

-4.1 

-3.8 

-0.2 

-2.8 

-2.9 

-3.4 

-2.0 

-2.5 

-0.9 

-3.6 

4.3 

7.6 

0.3 

-0.6 

-1.6 

0.1 

0.7 

-0.6 

1.5 

4.6 

3.3 

A32 

0.8 

2.0 

-2.5 

-5.7 

1.2 

-3.3 

3.0 

-3.1 

-8.8 

-5.0 

-4.1 

-5.4 

-2.2 

S2 

-3.1 

-4.6 

-0.2 

-1.0 

1.4 

2.7 

2.4 

2.6 

3.2 

3.5 

0.1 

-7.3 

-7.0 

A3 

-1.1 

-4.2 

2.5 

6.0 

-2.2 

-1.8 

1.7 

-3.9 

-4.9 

-2.5 

-8.0 

-5.8 

-3.4 

S3 

-3.9 

-5.7 

7.8 

6.7 

-3.5 

-4.5 

-3.1 

-0.3 

0.9 

1.4 

-7.2 

-11.3 

-10.3 

t Columns  G list  coefficients  of  the  velocity  anomalies  that  give  an  exact  match  to  geopotential  of  degree  2 and  3 
assuming  $v/£p  = — 6(knysec^(^cm^).  These  values  have  been  obtained  by  an  appropriate  perturbation  af  the  W* 
solutions  far  Regions  IIB,  III,  and  IV.  The  units  are  ny$ec. 
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where  is  the  density  kernel: 

Hii  - m t i hi + 3)  [(ri+i/a)i+3  - <Va)i+3]  ; (iv-5) 

and  are  the  observed  geopotential  coefficients  (for  C-P  we  have  used  the  value  of  Jef- 
freys1^ and  the  results  of  Gaposchkin11  for  the  remaining  9 coefficients)  and  a is  the  propor- 
tionality coefficient: 

<5v  = a6p 

Solutions  are  found  for  a wide  range  of  values  of  the  parameter  c*.  We  choose  the  value  of 
a that  requires  the  minimum  perturbation  in  the  observed  velocity  field. 

We  find  the  "best"  value  of  a to  be  -6  (km/sec)/(g/cm3),  assuming  that  the  density  per- 
turbations are  limited  to  Regions  11A,  111,  and  IV.  In  Fig.  1V-14,  we  compare  the  contours  of 
velocity  anomalies  obtained  from  the  original  W solution  (broken  lines)  with  the  solution  that 
gives  an  exact  match  to  the  gravity  field  for  a - ~6  (solid  lines).  The  spherical  harmonic  coef- 
ficients obtained  from  this  solution  are  incorporated  in  Table  1V-2,  columns  G.  The  changes 
in  the  pattern  of  velocity  anomalies  are  not  substantial.  It  should  be  remembered  that  the  ve- 
locity distribution  maps  contain  contributions  of  all  the  harmonics  with  i = 1,2,3.  Some  of  these 
harmonics  are  not  considered  in  comparison  with  the  gravity  field. 

The  negative  value  of  the  proportionality  coefficient  a suggests  several  possibilities  for  the 

origin  of  the  lower  mantle  velocity  anomalies:  (1)  sinking  of  dense,  low-velocity  material  into 

12  13 

the  lower  mantle  from  regions  of  lithospheric  subduction;  ' (2)  chemical  plumes  of  light, 

1 4 

high-velocity  material  originating  near  the  core-mantle  boundary;  (3)  temperature  differences 
and  perturbation  of  the  core-mantle  boundary  and  the  Earth's  surface  associated  with  mantle 
wide  convection,  and  (4)  static  chemical  heterogeneities  in  a non-convecting  mantle. 

The  first  three  alternatives,  all  involving  flow  in  the  lower  mantle,  may  be  complementary, 
but  act  on  a different  time  scale. 

There  appears  to  be  a westward  or  northwestward  translation  of  some  anomalies  with  re- 
spect to  the  pattern  obtained  for  the  innermost  shell.  In  particular,  the  direction  of  translation 
of  a large  negative  anomaly  under  the  Pacific  is  in  agreement  with  the  sense  of  motion  of  the 
Pacific  plate.  We  must  caution  the  reader,  however,  that  this  is  a highly  speculative  interpre- 
tation. If  correct,  it  would  favor  the  convective  hypotheses. 

A.  M.  Dziewonski 

B.  H.  Ilagert 

R.  J.  O'Connell t 

D.  INTERNAL  DISCONTINUITIES  IN  THE  EARTH'S  STRUCTURE 
AND  PERTURBATION  THEORY 

15-17 

Inversion  of  free  oscillation  data  has  led  to  the  currently  preferred  value  for  the  radius 

1 8 

of  the  core-mantle  boundary  (3485  ± 3 km).  This  value  has  since  been  strongly  supported  by 
1 9 

body  wave  studies.  This  "success"  maybe  in  retrospect  somewhat  surprising,  as  it  has  re- 
cently become  apparent  that  the  differential  kernels  used  to  predict  changes  in  eigenfrequencies 
of  spheroidal  and  radial  modes  due  to  perturbations  in  the  core  radius  were,  in  fact,  grossly 
erroneous. 


t Department  of  Geological  Sciences,  Harvard  University,  Cambridge,  MA  02138. 
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Backus  and  Gilbert^  were  the  first  to  propose  an  expression  for  perturbation  of  the  eigen- 
frequency  of  a normal  mode  due  to  a change  in  location  of  a discontinuity  (see  their  Appendix  B). 
Their  equation  (52),  hereafter  referred  to  as  L,  was  derived  by  what  is  now  known  to  be  an  in- 
complete application  of  the  Rayleigh  variational  principle.  The  L formula  was  criticized  by 
21 

Wiggins,  but  he  failed  to  provide  a rigorous  alternative  procedure.  In  Table  IV-3,  we  compare 
results  obtained  by  application  of  the  L formula  with  exact  calculations.  The  values  listed  cor- 
respond to  a relative  change  in  eigenfrequency  due  to  a 1-km  increase  in  the  radius  of  a discon- 
tinuity. The  "exact"  values  represent  the  relative  difference  between  the  eigenfrequencies  of 
the  perturbed  and  original  model  PEM-A  (Ref.  22)  computed  by  solving  the  appropriate  system 
of  differential  equations.  These  eigenfrequencies  are  determined  with  a relative  precision  of 


TABLE  IV-3 

RELATIVE  CHANGE  IN  EIGENFREQUENCY  (X  103)  DUE  TO  A 1-km  INCREASE 
IN  THE  RADIUS  OF  DISCONTINUITIES  IN  MODEL  PEM-A  (Ref. 22) 

Core-Mantle 

Boundory 

670-km  Disk 

Mode 

L 

Exact 

H 

L 

Exact 

H 

o 

V) 

o 

-0.0968 

-0.0955 

-0.0954 

-0.0092 

-0.0024 

-0.0024 

5S0 

0.0785 

-0.0781 

-0.0780 

-0.0254 

0.0158 

0.0157 

0S2 

-0.4652 

-0.2271 

-0.2270 

-0.0117 

0.0149 

0.0149 

0S 15 

-0.0012 

-0.0026 

-0.0026 

-0.1262 

0.0884 

0.0884 

0S30 

0.0000 

-0.0004 

-0.0005 

-0.0821 

0.1089 

0.1085 

0S60 

0.0000 

-0.0001 

-0.0001 

-0.0042 

0.0072 

0.0071 

1S8 

-0.6586 

-0.1077 

-0.1081 

-0.0187 

0.0162 

0.0162 

1S61 

0.0000 

-0.0001 

-0.0001 

-0.1260 

0.1276 

0.1271 

6S1 

0.0211 

-0.0133 

-0.0131 

-0.0317 

0.0241 

0.0240 

0T4 

-0.0345 

-0.0346 

-0.0345 

-0.0407 

0.0609 

0.0608 

0T15 

0.0000 

0.0000 

0.0000 

0.0251 

0.0940 

0.0938 

0T30 

0.0000 

0.0000 

0.0000 

0.0058 

0.0370 

0.0368 

0T60 

0.0000 

0.0000 

0.0000 

0.0003 

0.0032 

0.0032 

1 T43 

0.0000 

0.0000 

0.0000 

0.0819 

0.1376 

0.1371 

2T4 

0.1716 

0.1717 

0.1716 

-0.0776 

0.0697 

0.0694 

fColumns  L and  H refer  to  the  Lograngian  and  Homiltonion  representations,  respectively. 
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-7 

at  least  10  . The  corresponding  values  differ  significantly;  for  the  6 70 -km  discontinuity,  the 

L formula  yields  a wrong  sign  for  all  the  radial  and  spheroidal  modes  listed,  with  the  excep- 
tion of  qSq. 

This  result  was  communicated  to  several  scientists,  and  the  source  of  the  error  in  the  L 

23 

formula  was  detected  almost  immediately  by  Dr.  J.  H.  Woodhouse.  He  has  shown  that  a cor- 
rect application  of  the  Rayleigh  principle,  in  which  all  pertinent  variations  are  included,  leads 
to  an  expression  (formula  H)  that  is  significantly  different  from  the  L formula,  with  the  ex- 
ception of  those  boundaries  at  which  tractions  vanish  (free  surface,  core-mantle  boundary  for 

23  24 

toroidal  modes).  Algebraic  details  can  be  found  in  Woodhouse  and  Dziewonski  and  Sailor. 
Results  obtained  using  the  H formula  are  incorporated  in  Table  IV -3;  the  agreement  with 
"exact"  calculations  is  excellent. 

Note  that  for  the  toroidal  modes  ^T^  and  ^T^  both  L and  H formulas  give  identical  results 

at  the  core-mantle  boundary.  Also,  for  spheroidal  modes,  the  signs  are  the  same  even  though 

the  absolute  values  differ.  It  is  for  this  reason  that  inversion  of  a combined  set  of  over  1000 

1 7 

normal  modes  yielded,  after  several  iterations,  a correct  value  for  the  core  radius,  in  spite 
of  the  fact  that  the  incorrect  L formula  was  applied. 

Differential  kernels  for  perturbations  in  the  phase  velocity  due  to  a change  in  depth  to  a 
discontinuity  (such  as  Moho)  are  required  in  the  inversion  of  surface  wave  dispersion  data. 

For  an  increase  6h  in  the  depth  of  a discontinuity  located  in  the  original  model  at  z = — h,  the 
change  in  phase  velocity  c measured  at  a constant  frequency  o>  is 


(6c/  c) 


| £ 6h  [kK  + (iM  + pR]  + 


z=-h 


where  u is  the  group  velocity  and 
K = k2V2  - OzU)2  , 

M = i [(2 f>zU) 2 - (kV)2]  + (kU)2  - (SZV)2  + (kV)2  - (S  W)2  + (kW)2  , 

R = -w2(U2  + V2  + W2)  ; 

U and  V are  the  scalar  vertical  and  horizontal  displacement  functions  for  Rayleigh  waves,  WT 
is  the  displacement  function  for  Love  waves,  and  k is  the  wavenumber.  The  functions  U,  V, 
and  W are  normalized  such  that: 


■T 


p(U2  + V2)  dz  = 1 


-*  r 

_ oo 


pW  dz  = 1 


The  effect  of  changes  in  radii  of  discontinuities  is  also  important  in  evaluation  of  multiplet 

splitting  for  a laterally  heterogeneous  earth.  The  original  equations^ for  the  splitting  due 

to  ellipticity  were  derived  using  the  same  incorrect  approach  that  led  to  the  L formula.  Nu- 
2 7 

merical  calculations  based  on  this  approach  give  very  large  ellipticity  splitting  parameters 
for  earth  models  with  discontinuities  in  the  upper  mantle.  It  was  this  intuitively  unexpected 
result  that  originally  led  the  first  two  of  us  to  reexamination  of  the  entire  problem.  Dziewonski 
and  Sailor^  and  Dahlen^  have  now  derived  correct  equations  for  the  ellipticity  splitting.  Fur- 
ther theoretical  development  is  needed  for  the  general  case  of  a laterally  heterogeneous  earth 
model. 
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In  Figs.  IV -15  and  -16,  we  compare  for  three  earth  models  the  results  obtained  for  the  el- 

25 

lipticity  splitting  parameter  a defined  by  Dahlen  using  the  correct  (H)  and  incorrect  (L)  for- 
mulas for  contribution  from  discontinuities.  All  three  models  yield  numerical  results  which 
are  very  close  to  each  other  when  the  correct  equations  are  applied.  Differences  in  a in  the 
period  range  shown,  are  of  the  order  of  l/300  for  spheroidal  modes,  much  below  the  thickness 
of  the  line.  They  reach  1.5  percent  for  the  toroidal  modes  at  150  sec,  but  for  simplicity  we 
have  shown  only  a(T)  for  PEM-A  as  representative  of  all  three  models. 

Since  both  ellipticity  and  rotational  splitting  are,  for  practical  purposes,  essentially  inde- 
pendent of  the  details  of  an  acceptable  earth  model,  in  Table  IV -4  we  list  for  future  reference 


FUNCTIONS 

TABLE  IV-4 

Xi  AND  X2  (X103)  OF  DAHLEN27 
FOR  MODEL  PEM-A  (Ref. 22) 

COMPUTED 

Period 

(sec) 

Rayleigh  Wave 

Love  Wave 

co/uo 

X1 

X2 

c0/u0 

xi 

*2 

100 

1 .066 

-0.015 

0.617 

1.052 

0.014 

0.610 

125 

1.110 

-0.004 

0.656 

1.069 

0.023 

0.629 

150 

1.158 

0.004 

0.699 

1.086 

0.034 

0.650 

175 

1 .212 

0.007 

0.746 

1.103 

0.047 

0.672 

200 

1 .268 

0.004 

0.797 

1.120 

0.064 

0.695 

225 

1 .325 

-0.010 

0.846 

1.138 

0.083 

0.720 

250 

1.373 

-0.039 

0.889 

1.155 

0.102 

0.745 

275 

1 .405 

-0.089 

0.918 

1.171 

0.132 

0.770 

300 

1 .418 

-0.160 

0.930 

1.186 

0.161 

0.795 

325 

1 .414 

-0.249 

0.926 

1.199 

0.194 

0.819 

350 

1 .398 

-0.346 

0.912 

1 .211 

0.231 

0.842 

375 

1.378 

-0.443 

0.894 

1.221 

0.272 

0.865 

400 

1.357 

-0.531 

0.876 

1.229 

0.316 

0.886 

425 

1.338 

-0.605 

0.859 

1.235 

0.363 

0.906 

450 

1.320 

-0.662 

0.843 

1.239 

0.414 

0.925 

475 

1 .301 

-0.697 

0.826 

1 .241 

0.468 

0.942 

500 

1.279 

-0.705 

0.806 

1.242 

0.526 

0.959 

27 

functions  X^(T)  and  X2<T)  in  Eq.  (27)  of  Dahlen  for  the  apparent  length  of  a great  circular  path: 

L (T,  0)  = 2)ra  [l  - X,  (T)  cos  O - X,(T)  (1-3  cos2  9)] 
app  1 £ 

where  a is  the  average  earth  radius,  T is  the  period,  and  0 is  the  colatitude  of  the  pole  of  the 
great  circular  path  (0  = 0°  for  a wave  traveling  along  the  equator  from  west  to  east;  0 = 180° 
for  the  wave  traveling  in  the  opposite  direction),  and  X1  and  X^  are  functions  dependent  on  the 
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rotational  and  ellipticity  splitting,  respectively.  The  values, 
c0/V  have  been  computed  for  model  PEM-A  (Ref.  22). 


listed  together  with  the  ratios 


A.  M.  Dziewonski 
R.  V.  Sailort 
F.  A.  Dahlen^ 


E.  TRANSMISSION  HOLOGRAPHY  WITH  IMPULSIVE  SOURCES 

The  isolation  of  that  part  of  the  seismogram  due  to  scattering  provides  data  on  the  non- 
radially  symmetric  part  of  the  velocity  structure  and  an  accurate  picture  of  the  source  wave- 
form. In  the  high-frequency  limit,  the  path  effect  can  be  approximated  by  ray  theory  but  as 

much  as  50  percent  of  the  energy  in  a short -period  seismogram  cannot  be  so  explained  and  must 

29 

be  attributed  to  higher  order  scattering.  Chernov's  acoustic  scattering  model  has  been  ap- 

30 

plied  to  array  data  from  NORSAR,  LASA,  Garm,  and  California.  The  spatial  correlation 

2 

function  determined  for  LASA  and  NORSAR  reveals  that  the  scattering  is  not  due  to  isotropic 

random  inhomogeneities.  Additionally,  high-resolution  frequency-wavenumber  analysis  of  seis- 

3 1 

mic  energy  from  local  strip-mine  explosions  at  LASA  demonstrates  that  single  scattering  is 

not  the  cause  of  observed  near-field  seismic  signatures. 

An  alternative  to  statistical  studies  is  the  direct  inversion  of  travel-time  residuals  into 

32 

estimates  of  the  velocity  structure.  The  method  requires:  a large  number  of  events  whose 
ray  paths  uniformly  sample  the  structure  under  consideration,  an  initial  phase  assumption  for 
each  event,  and  a known  first-order  velocity  structure.  The  over-determined  inverse  problem 
is  numerically  difficult  to  solve  and  may  result  in  non-unique  solutions.  An  entirely  different 
approach  to  the  scattering  problem  is  to  downward  propagate  the  observed  wavefields.  This 
backward  propagation  technique  is  analogous  to  optical  holography  methods  and  can  be  expressed 
in  the  same  mathematical  form.  Application  of  the  holographic  technique  and  the  resulting  dif- 
ferences from  a ray  theory  inversion  can  be  illustrated  by  a homogeneous  half-space  model 
containing  a number  of  anomalous  velocity  structures.  In  Fig.  IV -17,  it  can  be  seen  that  the 
ray  path  A-A'  does  not  sample  the  anomalous  velocity  regions  1^  and  1^.  However,  the  recorded 
wavefield  will  contain  not  only  the  unperturbed  signal  but  also  the  scattered  wavefronts  from  1 
and  1^.  The  superposition  of  the  waves  from  all  the  sources,  both  virtual  and  primary,  can  be 

represented  by  a time-dependent  complex  amplitude  distribution  over  the  surface  plane  (x,  y). 

33 

The  amplitude  distribution  at  a given  time  can  be  represented  by: 

A(x,  y,  0)  = a(x,  y)  exp  [i<P(x,  y)  1 

The  relation  between  this  complex  amplitude  distribution  and  the  complex  amplitude  distribution 
B(x,  y,  Zq)  at  depth  z^  can  be  written  as: 

B(x,  y.  z0)  = T § \ AU  - ’J  • °)  ■ T(x,  y,  I , 77 ) d|  dT) 


where 


T = — exp  [-ikr]  cos  (n 


r) 


t Department  of  Geological  Sciences,  Harvard  University,  Cambridge,  MA  02138. 
t Department  of  Geological  and  Geophysical  Sciences,  Princeton  University,  Princeton,  NJ  08540. 
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n is  the  normal  vector  to  the  (x,  y)  plane,  and  r is  the  vector  connecting  the  point  (x,  y,  z^)  with 
the  point  ({  , 17 , 0) . 

It  can  be  seen  from  this  relation  that  each  point  (£  , v t 0)  maps  to  every  point  in  the  (x,  y,  zQ) 

plane  and  vice  versa.  This  is  essentially  Huygen's  principle  which  states  that  each  point  on  a 

plane  which  intersects  a wavefield  can  be  considered  a point  radiator.  Thus,  theoretically,  only 

one  recorded  wavefield  over  a plane  is  needed  to  recover  lateral  variations  in  velocity  beneath 

the  plane.  In  practice,  many  events  will  be  needed  to  average  the  effects  of  errors  in  the  data 

and  errors  due  to  undersampling  in  the  observation  plane. 

The  actual  mechanics  of  recovering  the  wavefield  at  various  depths  can  take  a number  of 

forms.  The  most  straightforward  technique  is  to  merely  replace  the  double  integral  in  the 

above  equation  and  evaluate  the  contributions  from  each  point  in  the  surface  plane.  This  brute- 

2 

force  technique  was  used  to  obtain  results  shown  earlier.  Another  promising  method  being 

34 

studied  utilizes  backward  wave  equation  propagation  based  on  the  theory  of  Claerbout  and 
35 

Landers  and  Claerbout.  For  these  techniques,  better  resolution  can  be  obtained  by  utilizing 
broadband  sources. 

In  conjunction  with  the  numerical  experiments,  a laboratory  experiment  using  impulsive 

sources  is  currently  nearing  completion.  The  apparatus  uses  a broadband  piezoelectric  source 

which  generates  2-MHz  pulses,  a movable  transducer  that  records  P-waves  and  epoxy  block 

velocity  models.  The  source  and  receiver  may  be  used  in  any  array  geometry  by  summing  the 

results  of  different  experiments.  The  laboratory  experiments  will  provide  data  about  wavefields 

produced  by  a number  of  known  inhomogeneities.  Experiments  involving  irregular  interfaces 

and  irregularly  shaped  velocity  anomalies  including  liquid-filled  cavities  similar  to  magma 

chambers  are  planned.  T , . 

r J.  bcheimer 

T.  E.  Landers 
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Fig.  IV-1.  Amplitude,  envelope,  and  exponential  traces  for  the  polarized  coda  of  event  1. 
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Fig.  IV -3.  Equal  area  projection  envelope 
peaks  of  2nd  minute  data. 


Fig.  IV-2.  Equal  area  projection  of  parti- 
cle motion  direction  for  peaks  of  envelope 
trace  for  the  1st  minute  data  of  event  1. 
The  rings  are  at  60°,  30°,  and  0°  emer- 
gence angle. 


Fig.  IV -4.  Equal  area  projection  envelope 
peaks  for  3rd  minute  data. 
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Fig.  IV-5.  Equal  area  projection  envelope 
peaks  for  4th  minute  data. 


Fig.  IV-6.  Equal  area  projection  for  4th 
minute  data  using  peaks  in  amplitude  trace. 
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Fig.  IV-8.  Radial  component  of  the  polarized  data  for  6 Novaya  Zemlya  events.  Plotting  gain 
is  the  same  as  for  Fig.  IV-7. 
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Fig.  IV-9.  P-to-P  scattering  region  for  phase  arriving  at  1 min.  17  sec  after 
initial  P phase.  Contours  indicate  epicenter  for  the  same  depth. 


73 


Tii-Mmil 


360 


270 


ISO 


Fig.  IV-10.  P-to-S  scattering  region  for  phase  arriving  at  2 min.  55  sec  after 
initial  P phase.  Contours  indicate  epicenter  for  the  same  depth. 
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Fig.  IV-11.  Mislocations  on  Scandinavian 
network  of  events  in  the  China-Russian 
border  region.  Head  and  tail  of  the  ar- 
rows represent  the  reference  and  the  ob- 
served location,  respectively.  The  ex- 
plosion site  in  Eastern  Kazakh  is  near 
(50  ° N,  78° E) . 


Fig.  IV-12.  Mislocations  on  Scandinavian  network  of  events  in  Central  America. 
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Fig.  IV-13.  Mislocations  on  Montana  network 
of  events  in  Bonin,  Mariana,  and  Ryukyu. 
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Fig.  IV-14.  Velocity  anomalies  in  (a)  Region  IIB  (1100  to  1500  km),  (b)  Region  III 
(1500  to  2200  km),  and  (c)  Region  IV  (2200  km  to  core-mantle  boundary).  Contours 
indicated  by  broken  line  represent  the  W solution  of  Table  IV-2;  solid  contours  cor- 
respond to  the  solution  perturbed  to  satisfy  the  gravity  data  (G).  The  results  in  this 
figure  correspond  to  inversion  for  a five-shell  model  of  693,495  data  corrected  for 
ellipticity  and  station  anomalies.  Contour  units  - m/ sec. 
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Fig.  IV-14.  Continued. 


77 


o 


& u 


j 1066  B / 


/ / 
/ / 

/ / 

' / 

/ 

/ 

"/ 


[-M31331 


\ 

\ 

\ 

\ \ 

\ \ 

\ \ 

\ N 
\ \ 

\ \ 

N \ 

\ N 
N S 
N 


1066  A , 1066  B ( PEM-A 


1066 A CONVENTIONAL  at  _ 


H ALGORITHM 

L ALGORITHM 


300  400 

PERIOO  (s*c) 


Fig.^IV-15.  Comparison  of  ellipticity  splitting 
parameters  for  the  fundamental  spheroidal 
mode  computed  for  three  earth  models.  Only 
the  algorithm  based  on  Hamiltonian  represen- 
tation (H  algorithm)  gives  correct  results. 
The  horizontal  line  represents  the  asymptotic 
value  (t  00 ) of  the  splitting  parameter. 


Fig.  IV-16.  Same  as  Fig.  IV-15,  but  for  the 
fundamental  toroidal  mode. 


Fig.  IV-17.  Model  of  a simple  medium 
illustrating  that  the  array  receives  in- 
formation from  inhomogeneities  that 
are  not  intercepted  by  the  geometrical 
ray  path.  The  inhomogeneities  are 
regions  of  higher  velocity  than  the  sur- 
rounding medium. 
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V.  MISCELLANEOUS  STUDIES 


A.  MODEL  ORDER  CRITERIA  APPLIED  TO  AUTOREGRESSIVE 
AND  NON-AUTOREGRESSIVE  PROCESSES 

Determining  the  proper  filter  length  to  use  is  crucial  to  estimating  the  spectrum  by  the 
maximum  entropy  method.  Several  estimating  procedures  for  determining  the  order  of  auto- 
regressive series  have  been  developed  and  the  equivalence  between  maximum  entropy  and 

1 

autoregressive  analysis  demonstrated.  Here  a study  of  these  procedures  for  both  autoregres- 
sive and  non-autoregressive  series  is  presented. 

Zero  mean  autoregressive  processes  have  the  following  form 


X.  = A . X,  . + A^X,  0 + . 
t 1 t - 1 Z t - Z 


. . + A X,  + e. 
p t-p  t 


where  et  is  zero  mean  and  uncorrelated,  is  referred  to  as  the  innovation.  Thus,  the  time 

series  X is  predictable  from  a linear  combination  of  its  pass  values  to  a random  error.  Var- 
1 2 

ious  ways  of  estimating  the  A.'s  are  available  and  will  not  be  discussed  here.  The  determina- 
tion of  the  value  of  the  order  of  the  process,  p,  is  a major  problem  in  the  analysis  and  is  the 
subject  of  this  research. 

3 

The  following  methods  have  been  proposed  to  estimate  p.  Akaike  suggested  that  the  min- 
imum of  the  average  error  due  to  estimating  the  autoregressive  coefficients  and  the  innovation 
for  one  step  prediction  gives  the  model  order,  p.  The  criterion,  called  the  final  prediction 
error,  FPE,  to  be  minimized  is 


FPE(M)  = P [N  + (M  + 1 ) ]/[ N - (M  + 1)] 


where  P^  is  the  residual  squared  error  for  an  Mth  length  filter  and  N is  the  length  of  X.  The 
M for  which  FPE  is  a minimum  is  taken  as  p.  Alternatively,  Akaike4  suggested  minimizing 
the  log  likelihood  of  the  innovation  variance  as  a function  of  filter  length  to  find  p.  This  crite- 
rion, called  the  information  theoretic  criterion,  AIC,  is  estimated  by 

AIC(M)  = In  (PM)  + 2M/N  . 


Again,  the  M for  which  AIC  is  minimized  is  taken  as  p.  The  third  method  considered  here 
was  proposed  by  Parzen3  and  is  known  as  the  autoregressive  transfer  function  criterion,  CAT. 
The  order  p is  given  where  the  estimate  of  the  difference  of  the  mean  square  error  between 
the  true  filter,  which  exactly  gives  the  innovation,  and  estimated  filter  is  a minimum.  Parzen 
showed  that  this  difference  can  be  calculated,  without  explicitly  knowing  the  exact  infinite  filter, 
by 


CAT(M) 


whe  re 

P.  = N/(N  - j)  P . 

The  following  examples  were  generated  to  compare  the  method  on  both  auto-  and  non- 

1 

autoregressive  series.  For  each  example  the  Burg  method  was  used  to  determine  the  A.'s. 
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For  plotting  purposes,  noting  that  log  FPE  asymptotically  approaches  AIC,  we  define 
F(M)  = log10  {log10  lFPE(M) ] - log1Q  f FPE(p) ] + 1} 

A(M)  = log10  [AIC(M)  - AIC(p)  + 1) 

C(M)  = l°g10  (CAT(M)  - CAT(p)  + 1] 

such  that  the  value  of  each  criterion  at  the  predicted  order  number,  p,  is  0. 

The  first  time  series  considered  is  one  cycle  of  the  complex  harmonic  exp[v^~ 1 (27r/l00)  t] 
shown  in  part  (a)  of  Fig.  V-l.  Part  (b)  of  the  same  figure  shows  the  three  order  criteria,  FPE, 
AIC,  and  CAT.  Each  curve  has  its  minimum  at  the  same  value,  the  CAT  being  much  more 
robust.  Using  that  value  for  p results  in  the  correct  spectrum  shown  in  part  (c).  Parts  (d) 
and  (e)  show  the  order  criteria  and  spectra,  respectively,  for  the  real  part  of  the  time  series 
and  parts  (f)  and  (g),  the  same  for  the  imaginary  part.  In  each  case  the  same  order  number  is 
indicated  and  the  corresponding  correct  spectrum  generated.  Again,  the  CAT  criterion  is  more 
robust  even  though  the  behavior  is  somewhat  artificial  due  to  the  use  of  logarithms. 

Figure  V-2  shows  the  results  for  the  same  time  series  plus  50  percent  uniformly  distrib- 
uted white  noise  added  independently  to  the  real  and  imaginary  parts.  Part  (a)  is  the  time 
series,  part  (b)  the  order  criteria  for  the  complex  case,  and  part  (c)  the  spectrum  assuming 
the  predicted  order  number.  Parts  (d),  (e),  and  (f)  are  the  order  criteria  and  spectra,  based 
on  CAT  and  AIC  and  the  FPE  minimas,  respectively,  for  the  real  part  of  the  time  series.  The 
CAT  and  AIC  gave  the  same  lower  value  for  p but  failed  to  produce  the  correct  spectrum,  while 
the  FPE  criterion  gave  the  higher  value  of  p and  the  correct  spectrum.  The  FPE  does  have  a 
local  minimum  of  almost  the  same  value  as  the  absolute  minimum  where  the  AIC  and  CAT 
curves  had  a minimum.  In  parts  (g),  (h),  (i),  and  (j),  the  order  criterion  and  spectra  at  p, 

2p  = p (at  the  FPE  real  value),  and  3p  for  the  imaginary  part  of  the  time  series  are  shown, 
respectively.  In  this  instance,  the  predicted  value  of  p was  the  same  for  the  AIC,  FPE,  and 
CAT  criteria;  however,  that  value  gave  a poor  spectral  estimate.  The  spectra  for  a filter 
length  of  3p  gave  the  best  result. 

The  second  experiment  used  the  Ricker  wavelet 

ft2  2 

R(t)  = — exp  [-7r(ktr) 

whose  power  spectrum  is 

R(f)  = <^-)4  exp  [-27r(^)2] 

as  the  time  series  whose  order  was  to  be  estimated  and  corresponding  spectrum  computed.  This 
function  is  not  autoregressive.  The  top  curves  of  Fig.  V-3  show  the  three  symmetric  portions 
of  the  wavelet  for  k = 1,  2,  and  10,  respectively.  The  next  three  curves  show  the  filter  length 
criteria  for  the  corresponding  wavelet  above.  For  each  wavelet  the  AIC,  FPE,  and  CAT  had 
the  same  value,  although  it  increased  as  the  wavelet  became  compacted.  The  next  set  of  curves 
are  the  exact  log  and  linear  spectra  and  the  last  set  of  curves  are  the  log  and  linear  spectra 
computed  using  the  predicted  filter  length  criteria.  In  each  case  acceptable  spectra  were 
computed. 
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To  summarize  the  results,  we  note  that  in  the  case  where  little  noise  is  present,  AIC,  FPE, 
and  CAT  give  an  order  number  that  produces  an  acceptable  spectra  for  both  an  auto-  and  a non- 
autoregressive  time  series.  When  the  innovation  was  large,  the  criteria  gave  inconsistent  re- 
sults and  for  the  case  we  examined,  larger  filter  lengths  produced  the  correct  spectra. 

It  should  be  noted  that  the  criteria  for  order  number  are  based  on  minimizing  or  predicting 
the  innovation  and,  thus,  they  produce  not  what  one  usually  wants  to  observe  (the  spectra  with- 
out the  influence  of  any  noise  that  is  present)  but  a spectra  including  the  effects  of  noise.  In 
addition,  they  are  influenced  by  the  technique  used  to  get  the  autoregressive  coefficients.  We 
note  that  using  the  Burg  technique,  we  obtained  a large  filter  length  when  exp  1*1^1  Qt ] can  be 
exactly  predicted  by  the  one  point  prediction  error  wavelet  {l,  — expl*f-l  (Jq]}.  Work  is  con- 
tinuing to  generalize  the  observations  made  in  this  report. 

T.  E.  Landers 


B.  APPLICATION  OF  THE  ENERGY-MOMENT  TENSOR  RELATION 

TO  FINITE  DURATION  SOURCES 

Previously,^  we  derived  an  expression  relating  the  total  radiated  seismic  energy  from  an 
arbitrary  point  source  mechanism  to  the  moment  tensor  of  the  source.  This  expression  is 
formulated  in  terms  of  the  free  oscillation  normal  modes  of  the  earth  which  limits  its  applica- 
ble bandwidth  to  that  of  a given  normal  mode  catalog.  We  showed  results  for  two  large  earth- 
quakes (Chile,  1960;  Alaska,  1964)  for  each  of  two  source  time  functions:  a step  and  a ramp 
increase  in  applied  moment.  In  this  report,  we  show  results  from  a third  source  time  function, 
a half- sine -wave  increase  in  applied  moment,  and  discuss  the  physical  implications  and  limita- 
tions on  the  radiated  energy. 

The  three  source  time  functions  are  shown  in  Fig.  V-4.  The  first  two  are  as  previously 

7 

described;  the  third  is  the  same  as  that  used  by  Kanamori  and  Anderson.  From  these  functions, 
energy-moment  tensor  relations  may  be  easily  derived  for  the  total  radiated  energy.  These  are: 
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where  MQ  is  the  final  value  of  the  moment  tensor,  Sn(rg)  is  the  strain  tensor  of  the  n**1  normal 

mode  evaluated  at  the  source  position,  F is  the  normalization  integral  (previously  described), 

th  ^ 

w is  the  angular  frequency  of  the  n n normal  mode,  and  T is  the  source  duration, 
n 8 

We  present  below  results  for  the  I960  Chilean  earthquake  based  on  Plafker  and  Savage's 

29 

fault  plane  solution  with  a seismic  moment  of  8.6  X 10  dyn-cm.  These  are  computed  from  a 

9 

normal  mode  catalog  originally  supplied  to  us  by  Buland  and  Gilbert  and  now  extended  to  pro- 
vide complete  coverage  down  to  a 45 -sec  period  for  both  poloidal  and  toroidal  modes. 
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Figures  V-5  to  -7  show  the  effect  of  changing  the  source  time  function  on  the  spectral  dis- 
tribution of  energy  in  the  toroidal  modes.  The  results  are  similar  for  the  poloidal  modes.  These 
plots  are  cumulative  energy  in  0.004  rad/sec  intervals  of  angular  frequency.  The  period  range 
of  the  plots  is  from  18,000  sec  at  the  left-hand  side  to  45  sec  at  the  right-hand  side.  This  was 
chosen  to  encompass  the  limits  of  our  normal  mode  catalog. 

The  results  in  these  figures  show  that  the  step  function  source  produces  an  energy  spectrum 
which  continues  to  increase  as  cj1,85  at  short  periods.  This  is  clearly  impossible  because  it 
implies  an  infinite  radiated  energy.  The  ramp  function  source  results  decrease  as  w at 

short  periods  which,  although  they  are  certainly  preferable  to  those  of  the  step  function  source, 
still  fail  to  guarantee  a finite  radiated  energy.  Only  the  third  source  time  function,  a half- sine  - 
wave  increase,  produces  an  energy  spectrum  which  is  convergent.  The  drop-off  in  radiated 
energy  varies  as  w for  this  source.  We  must  therefore  reject  the  step  and  ramp  source 
time  functions  as  physically  unrealizable  because  of  their  asymptotic  behavior  at  short  periods. 

Because  of  the  difficulties  of  measuring  the  amplitude  decay  of  spectra  from  their  enve- 
lopes, particularly  when  the  bandwidth  is  as  limited  as  ours,  we  cannot  give  any  certainty  to 
the  values  of  the  exponents  given  above.  Experience  suggests  the  fractional  error  in  the  slope 
of  a line  measured  off  a graph  to  be  bounded  by  twice  the  largest  fractional  error  in  the  ordi- 
nates. If  this  is  true,  then  we  expect  the  values  of  the  exponents  to  be  accurate  to  within  10  per- 
cent. Thus,  we  cannot  say  for  sure  whether  the  exponent  for  the  step  function  frequency  depen- 
dence, 1.85,  precludes  a frequency  squared  behavior.  On  the  other  hand,  the  half-sine -wave 
source,  dropping  off  in  frequency  with  a measured  exponent  of  2.6  seems  to  rule  out  either 
frequency  squared  or  frequency  cubed  behavior. 

The  dependence  of  the  radiated  energy  on  the  source  duration  for  the  half- sine -wave  source 
time  function  is  shown  in  Fig.  V-8.  As  can  be  seen,  there  is  a dramatic  decrease  in  the  total 
radiated  energy  on  the  order  of  a factor  of  400  between  0-  and  200-sec  source  duration.  Further- 
more, this  decrease  is  logarithmically  nonlinear;  dropping  off  most  rapidly  at  short  source 
durations.  The  relative  distribution  of  radiated  energy  into  poloidal  and  toroidal  components  is 
also  dependent  on  source  duration.  v At  short  durations,  the  poloidal  component  dominates;  at 
long  durations,  the  radiated  energy  approaches  equipartition  between  poloidal  and  toroidal 
components. 

It  is  apparent  from  the  energy-moment  tensor  relation  that  the  partition  of  energy  between 
the  toroidals  and  poloidals  will  depend  on  the  fault  plane  solution  and  the  source  depth.  Since 
we  have  studied^  only  large  thrust  events  at  shallow  depths  (i.e.,  Alaska,  1964  and  Chile,  I960), 
our  results  cannot  be  extended  to  deep  events  of  different  mechanisms  without  modification. 

We  feel  the  physical  significance  of  this  phenomenon  lies  in  the  nature  of  the  forces  resist- 
ing the  earthquake  faulting  process.  Short  source  duration  (i.e.,  rapid  release  of  tectonic 
stress)  implies  there  are  low  frictional  or  deformational  forces  resisting  the  rupture  formation. 
Thus,  relatively  little  energy  is  expended  in  overcoming  those  resisting  forces.  Long  source 
duration  implies  the  opposite:  that  the  resisting  forces  are  large  and  therefore  consume  much 
of  the  tectonic  strain  energy  available  to  produce  the  rupture.  In  the  former  case,  the  seismic 
efficiency  is  high.  In  the  latter  case,  it  is  low  because  the  bulk  of  the  energy  is  wasted  in  non- 
seismic  processes.  The  results  of  Fig.  V-8  show  the  extent  of  this  variation  which,  for  our 
model,  spans  almost  2 surface  wave  magnitude  units. 

D.  W.  McCowan 
A.  M.  Dziewonski 
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C.  CONTRIBUTION  OF  LARGE  EARTHQUAKES  TO  EXCITATION  OF  POLAR 
MOTION  FOR  YEARS  1901  TO  1970 

Since  its  discovery  in  1891,  the  Eulerian  free  nutation  of  the  earth  (properly  called  the 
Chandler  wobble)  has  been  of  considerable  interest  to  astronomers,  meteorologists,  oceanog- 
raphers, and  solid -earth  geophysicists.  Throughout  this  time,  the  mechanisms  exciting  the 
wobble  have  remained  elusive,  and  although  atmospheric  motions,  core-mantle  interactions 
and  earthquakes  have  all  been  considered,  none  of  these  has  been  either  clearly  identified  as 
the  principal  source  of  excitation  or  entirely  rejected  as  a partial  contributor. 

In  this  report,  we  consider  the  cumulative  effect  of  major  earthquakes  that  occurred  be- 
tween 1901  and  1970  as  a source  of  excitation  of  the  wobble.  Previous  investigations  of  the 

role  of  earthquakes  in  exciting  the  wobble  have  yielded  no  clear  conclusions.  Although  Munk 
1 0 

and  MacDonald,  among  others,  had  concluded  that  displacements  associated  with  earthquakes 

1 1 

were  too  small  to  excite  the  wobble,  Smylie  and  Mansinha  presented  evidence  for  a correla- 
tion between  times  of  occurrence  of  major  earthquakes  and  changes  in  the  direction  of  polar 
motion.  This  suggested  that  earthquakes  could  be  responsible  for  a sizable  fraction  of  the  ex- 
citation and  revived  interest  in  the  subject  which  has  become  one  of  the  more  controversial 

topics  in  the  geophysical  literature. 

12  13 

Subsequent  studies  * concentrated  on  the  correlation  between  individual  earthquakes  and 
short-term  changes  in  polar  motion.  Because  of  the  considerable  uncertainty  in  the  exact  loca- 
tion of  the  pole  at  any  time,  a definite  correlation  between  polar  motion  and  earthquakes  could 

1415  16 

not  be  firmly  established.  ' In  addition,  Dahlen  concluded  that  earthquakes  could  not  ac- 
count for  the  observed  amplitude  of  the  wobble. 

Nevertheless,  there  is  still  highly  suggestive  evidence  of  a correlation  between  global 

17-19 

seismicity  and  changes  in  the  amplitude  of  the  wobble  over  the  period  of  the  last  70  years. 

This  suggests  that  the  cumulative  effect  of  earthquakes  may  account  for  long-term  variations 
in  the  amplitudes  of  the  wobble,  even  though  the  effect  of  any  individual  earthquake  is  too  small 
to  be  detected  above  the  noise  in  the  polar  motion  data.  In  order  to  investigate  this  possibility, 
we  calculate  the  polar  motion  expected  from  major  earthquakes  over  the  last  70  years  and  com- 
pare it  with  actual  polar  motion  data  over  the  same  period.  This  requires  specification  of  the 
time,  location,  and  source  parameters  of  the  earthquakes,  as  well  as  the  calculation  of  the 

change  of  the  earth’s  inertia  tensor  caused  by  an  earthquake. 

20 

In  our  previous  report,  we  have  described  a method  of  calculation  of  the  polar  shift  from 
an  earthquake  of  known  location,  source  mechanism,  and  seismic  moment.  We  also  have  dem- 
onstrated how  the  knowledge  of  relative  motions  of  plates  can  be  used  to  infer  the  geometry  of 
a fault  plane  solution  of  an  earthquake  occurring  at  a particular  type  of  the  boundary  between 
two  plates. 

The  234  earthquakes  with  magnitudes  M >7.8  used  to  calculate  the  excitation  of  the 

S 21 

Chandler  wobble  are  taken  from  the  catalog  compiled  by  Duda  and  supplemented  by  informa- 
tion from  the  Bulletin  of  the  Seismological  Society  of  America  for  the  years  1965-1970.  Since 
it  is  the  seismic  moment  MQ  that  determines  the  static  displacement,  we  have  used  an  approx- 
imate empirical  relation  to  determine  the  moments  from  the  reported  magnitudes.  Figure  V-9, 
after  Chinnery  and  North,  shows  the  relation  between  MQ  and  Mg  for  a number  of  earthquakes 
for  which  both  have  been  independently  determined.  The  solid  line  is  the  relation  between  mo- 
ment and  magnitude  proposed  by  Chinnery  and  North.  The  dashed  line  shows  the  relation  we 
have  used  to  determine  the  moments;  it  obviously  is  a reasonable  approximation  to  the  data  for 
magnitudes  greater  than  7.5.  The  remaining  source  parameters  were  determined  from  the 
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TABLE  V-l 

TIMES  AND  LOCATIONS  OF  THE  30  LARGEST  EARTHQUAKES 
CONTRIBUTING  TO  THE  EXCITATION  OF  THE  CHANDLER  WOBBLft 

No. 

Year 

Month 

Day 

Latitude 

Longitude 

Depth 

(km) 

Magnitude 

Shift 

(0.01  in.) 

Direction 

(°E) 

ALOD 

(psec) 

1 

1901 

Aug. 

9 

-22.0 

170.0 

s 

8.4 

0.60 

336.2 

-7.  02 

2 

1902 

Sept. 

23 

16.0 

-93.0 

s 

8.4 

0. 49 

296.3 

-5.83 

3 

1905 

Apr. 

4 

33.0 

76.0 

s 

8.6 

1.62 

89.4 

-2.60 

4 

1906 

Jan. 

21 

34.0 

138.0 

340 

8.4 

0.  83 

101.4 

1.84 

5 

1906 

Jan. 

31 

1.0 

-81.5 

s 

8.9 

1.39 

218.0 

-160. 17 

6 

1906 

Aug. 

17 

-33.0 

-72.0 

s 

8.6 

2.56 

108.0 

-16.09 

7 

1906 

Sept. 

14 

-7.0 

149.0 

s 

8.4 

0.79 

239.0 

-0.03 

8 

1910 

June 

16 

-19.0 

169.5 

100 

8.6 

1.78 

340.7 

-26.53 

9 

1911 

June 

15 

29.0 

129.0 

160 

8.7 

4.86 

105.7 

-17.96 

10 

1914 

Nav. 

24 

22.0 

143.0 

no 

8.7 

3.07 

160.6 

-43.  10 

11 

1917 

May 

1 

-29.0 

-177.0 

60 

8.6 

2.33 

9.0 

-18. 26 

12 

1917 

June 

26 

-15.5 

-173.0 

s 

8.7 

4.33 

277.3 

-0.69 

13 

1919 

Apr. 

30 

-19.0 

-172.5 

s 

8.4 

0.55 

19.6 

-8.  13 

14 

1922 

Nav. 

11 

-28.5 

-70.0 

s 

8.4 

0.  72 

110.8 

-6.  15 

15 

1923 

Feb. 

3 

54.0 

161.0 

s 

8.4 

0.71 

146.  1 

2.07 

16 

1929 

Mar. 

7 

51.0 

- 1 70.  0 

60 

8.6 

2.  10 

171.9 

6.95 

17 

1934 

Jan. 

15 

26.5 

86.5 

s 

8.4 

0.53 

95.3 

-3.40 

18 

1938 

Nov. 

10 

55.5 

-158.0 

s 

8.7 

3.  89 

189.  1 

17.  54 

19 

1939 

Dec. 

21 

0.0 

123.0 

150 

8.6 

0.  84 

106.7 

-4.36 

20 

1940 

May 

24 

-10.5 

-77.0 

60 

8.4 

0.33 

73.9 

-8.  76 

21 

1941 

June 

26 

12.5 

92.5 

60 

8.7 

4.30 

178.4 

-13.57 

22 

1941 

Nav. 

25 

37.5 

-18.5 

s 

8.4 

0.53 

323.6 

-0.27 

23 

1942 

Aug. 

24 

-15.0 

-76.0 

60 

8.6 

1.30 

75.6 

-22.71 

24 

1946 

Dec. 

20 

32.5 

134.5 

s 

8.4 

0.  77 

116.8 

-3.20 

25 

1950 

Aug. 

15 

28.5 

96.5 

s 

8.7 

2.86 

122.3 

-11.83 

26 

1952 

Mar. 

4 

42.5 

143.0 

s 

8.6 

2.  77 

131.2 

-5.00 

27 

1952 

Nav. 

4 

52.8 

159.5 

s 

8.4 

0.  82 

147.2 

1.69 

28 

1958 

Nav. 

6 

44.5 

148.5 

75 

8.7 

4.50 

130.4 

-2.47 

29 

1960 

May 

22 

-39.5 

-74.5 

s 

8.6 

2. 56 

109.2 

-9.80 

30 

1964 

Mar. 

28 

61.1 

-147.6 

20 

8.5 

1.11 

202.7 

7. 92 

|The  magnitude  and  direction  af  the  pole  shift  as  well  as  the  change  in  length  af  day  from  each  earthquake 
is  shown.  Depths  designed  s were  assigned  a standard  value  af  33  km.  Earth  model  PEM-A  (Ref.  24)  was 
used  to  calculate  the  pole  shifts. 
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plate  tectonic  model  using  the  data  of  Solomon  et  aL  for  the  positions  of  the  boundaries  and 
relative  motions  of  plates. 

The  calculated  pole  shift  (in  0". 01)  and  direction  (degrees  east  of  Greenwich)  are  given  in 
Table  V-l  for  the  30  largest  earthquakes  of  those  considered.  Also  shown  is  the  change  in 
length  of  day  (ALOD,  in  microseconds)  resulting  from  the  change  in  the  polar  moment  of  inertia 
of  the  earth. 

From  the  calculated  pole  shifts,  we  can  compute  a synthetic  polar  wobble  curve  that  would 
result  from  this  source  of  excitation  alone.  For  this,  we  also  require  the  starting  position 
m(t^)  of  the  pole  in  January  1901  and  the  frequency  o>q  and  damping  Q of  the  Chandler  wobble. 
The  synthetic  wobble  m(t)*is  computed  from- 

m(t)  = m(t0)  exp  [iw(t  - tQ)]  + £ skH(t  -tk)-  {1-(1+ir)  exp 

k 

where  u>  is  the  complex  circular  frequency  of  the  Chandler  wobble:  u>  = cjq(  1 + i/2Q),  s^  are 

the  complex  pole  shifts  due  to  earthquakes  at  times  t^,  ^ is  the  frequency  of  the  diurnal  rota- 
tion and  H is  the  unit  step  function.  The  complex  function  m(t)  is  referred  to  the  coordinate 
system  in  which  the  real  axis  corresponds  to  the  Greenwich  meridian  and  the  imaginary  axis 
points  toward  90 °E. 

Figure  V-10  shows  the  real  part  (m  ) of  three  such  synthetic  wobble  curves  obtained  using 

cJq  = 5.287  rad/yr  after  Wilson^  and  three  different  values  of  Q.  Shown  for  comparison  is  a 

smoothed  version  of  the  International  Latitude  Service  (ILS)  polar  motion  data,  with  the  annual 

25 

term  removed,  also  from  Wilson.  The  initial  pole  location  is  the  same  for  all  curves,  and 
each  curve  has  been  bandpass  filtered  from  0.3  to  0.98  cycles/yr  to  facilitate  comparison.  The 
ticks  on  the  top  of  the  figure  mark  the  times  of  the  earthquakes  in  Table  V-l. 

The  general  trend  of  the  amplitude  of  the  wobble  can  be  seen  in  the  ILS  data.  The  observed 
amplitude  was  large  in  1910,  decreased  fairly  uniformly  from  1910  to  1920,  and  remained  small 
until  1940  when  it  began  to  increase  again  to  reach  the  maximum  value  by  1954.  Since  then, 
there  has  been  an  irregular  decline  in  amplitude.  The  synthetic  wobble  time  series,  particu- 
larly for  Q = 100,  reproduce  the  principal  features  of  these  amplitude  changes.  The  decrease  in 
the  synthetic  wobble  is  due  to  the  cumulative  effect  of  several  large  earthquakes  from  1914  to 
1929.  This  occurs  even  for  the  case  with  Q = 10,000  and  it  is  thus  apparent  that  earthquakes 
can  reduce  the  amplitude  of  the  wobble  and  a rapid  decay  of  amplitude  need  not  be  attributed  to 
intrinsic  dissipation.  Similarly,  the  increase  in  the  synthetic  wobble  amplitude  is  due  to  large 
earthquakes  from  1938  through  195  0 which  gradually  re -excite  the  wobble.  We  should  note  that 
the  decline  in  amplitude  observed  since  1954  is  not  reproduced  in  the  synthetic  wobble  curves; 
this  may  be  an  effect  of  other  sources  of  excitation. 

The  results  of  this  study  indicate  that  the  cumulative  effects  of  major  earthquakes  can 
account  for  at  least  a large  part  of  the  excitation  of  the  Chandler  wobble  and,  in  particular,  for 

the  long-term  variations  of  the  amplitude  of  the  wobble.  This,  together  with  recent  estimates 

25 

of  the  atmospheric  excitation  of  the  wobble,  indicates  that  these  two  mechanisms  may  well 
account  for  the  entire  excitation.  If  this  is  verified  by  more  detailed  analysis,  it  should  per- 
mit more  accurate  determination  of  the  period  and  damping  of  the  wobble,  and  of  the  effective 
elastic  properties  of  the  earth  at  very  long  periods. 

R.  J.  O' Connell  t 
A.  M.  Dziewonski 

t Department  of  Geological  Sciences,  Harvard  University,  Cambridge,  MA  02138. 
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Fig.  V-l.  Error  analysis  for  one  cycle  of  a complex  harmonic:  (a)  Complex  harmonic 

exp  [\f^(27r/l00)t ];  (b)  AIC,  CAT,  and  FPE  for  (a);  (c)  spectrum  at  minimum  in  (b), 
(d)  AIC,  CAT,  and  FPE  for  real  part  of  (a);  (e)  spectrum  at  minimum  in  (d);  (f)  AIC, 
CAT,  and  FPE  for  imaginary  part  of  (a),  and  (g)  spectrum  at  minimum  in  (f). 
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Fig.  V-2.  Error  analysis  for  one  cycle  of  a complex  harmonic  with  white  noise:  (a)  Complex 
harmonic  exp [N/^l(27r/l00)t]  plus  50  percent  noise;  (b)  AIC,  CAT,  and  FFE  for  (a);  (c)  spec- 
trum  for  minimum  in  (b);  (d)  AIC,  CAT,  and  FPE  for  real  part  of  (a);  (e)  spectrum  for  min- 
imum given  by  AIC  and  CAT  in  (d);  (f)  spectrum  for  minimum  given  by  FPE  in  (d);  (g)  AIC, 
CAT,  and  FPE  for  imaginary  part  of  (a);  (h)  spectrum  for  minimum  in  (g);  (i)  spectrum  for 
twice  minimum  in  (g);  and  (j)  spectrum  for  three  times  minimum  in  (g). 
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Fig.  V-3.  Error  analysis  for  the  non-autoregressive  Ricker  wavelet:  (a)  Ricker  wavelet  k = 1;  (b)  AIC, 

FPE,  and  CAT  for  (a);  (c)  exact  spectrum  of  (a);  (d)  spectrum  of  (a)  based  on  minimum  in  (b);  (e),  (f), 

(g),  and  (h)  same  as  (a-d)  for  k = 2;  and  (i),  (j),  (k),  and  (1)  same  as  (a-d)  for  k = 10. 
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Fig.  V-4.  Three  different  source  time  functions. 
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Fig.  V-5.  Spectral  distribution  of  radiated 
energy  for  a step  source  time  function  of 
applied  moment.  Results  for  the  I960  Chil- 
ean earthquake  using  Plafker  and  Savage's 
(1970)  fault-plane  solution.  180-sec  source 
at  79.4  km  depth. 
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Fig.  V-6.  Spectral  distribution  of  radiated 
energy  for  a ramp  source  time  function  of 
applied  moment.  Results  for  the  I960  Chil- 
ean earthquake  using  Plafker  and  Savage’s 
(1970)  fault-plane  solution.  180-sec  source 
at  7 9.4  km  depth. 


Fig.  V-7.  Spectral  distribution  of  radiated 
energy  for  a half- sine -wave  source  time 
function  of  applied  moment.  Results  for  the 
I960  Chilean  earthquake  using  Plafker  and 
Savage’s  (197  0)  fault-plane  solution.  180-sec 
source  at  79.4  km  depth. 
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Fig.  V-8.  Energy  as  a function  of  source 
duration  for  the  half- sine -wave  source  time 
function.  Results  for  the  1970  Chilean  earth- 
quake using  Plafker  and  Savage's  (1970)  fault- 
plane  solution.  Source  at  79.4  km  depth. 
Magnitude  from  Gutenberg-Richter  relation. 
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Fig.  V-9.  Relation  between  seismic  moment  Mo  and  magnitude  Ms 
after  Chinnery  and  North.22  The  dashed  line  gives  the  relation  used 
to  calculate  the  pole  shifts  from  earthquakes. 
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Fig.  V-10.  Simulated  polar  motion  curves  resulting  from  excitation  by  large 
earthquakes;  the  damping  for  each  is  determined  by  Qs.  ILS  is  smoothed 
observed  polar  motion  with  anneal  term  removed.24  Only  the  component  along 
the  Greenwich  meridian  is  shown. 
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VI.  DATA  SYSTEMS 


A.  SEISMIC  DATA  AND  INFORMATION  SYSTEM 

A seismic  information  and  data  system  has  been  under  design  and  development  for  the  last 
year  and  a half.  The  progress  in  implementing  the  system  has  been  described  in  several  of  the 
more  recent  Semiannual  Technical  Summaries  issued  during  that  period.  This  present  report 
will  be  a project  milestone  in  that  it  will  mark  both  the  end  of  the  development  phase  of  the  sys- 
tem and  the  beginning  of  the  information  distribution  phase.  A complete  report  describing  the 
project  is  in  preparation  and  should  be  ready  for  distribution  during  the  third  quarter  of  1976. 
What  is  presented  in  this  short  report  is  a general  overview  of  the  system  and  the  various  files 
that  make  up  the  directory  (IWWSS)  at  ISIC  on  the  ARPANET. 

There  are  two  reasons  for  developing  such  a seismic  data  and  information  system,  and 
both  reasons  are  equally  important  for  the  development  of  this  data  base.  One  reason  was  the 
complexity  of  the  seismic  research  observatory  program  combined  with  the  previously  installed 
large  seismic  arrays  and  the  high-gain  long-period  installations.  These  seismic  systems  com- 
bined to  produce  a highly  sophisticated  system  of  digital  recording  seismic  stations.  This  sys- 
tem needed  a logically  structured  description  that  could  be  made  available  to  users  of  the  data 
generated  by  this  seismic  network.  A second  reason  was  the  development  of  a sophisticated 
text  editing  and  management  system  called  NLS  (oN-Eine  System)  developed  by  the  Augmentation 
Research  Center  of  Stanford  Research  Institute.  This  system  was  made  available  over  the 
ARPANET  to  a very  diverse  community  of  users.  NLS  was  to  provide  the  basic  tool  necessary 
to  build  the  data  base. 

The  Seismic  Data  and  Information  System  (SDIS)  is  an  information  data  base  that  describes 
in  detail  the  Seismic  Research  Observatory  (SRO)  program  development  and  the  high-gain  long- 
period  program.  The  older  large  seismic  arrays  are  also  covered;  however,  not  in  such  detail 
as  are  the  newer  components  of  the  seismic  network.  Additional  information  is  contained  in  the 
SDIS  describing  the  mass  storage  device,  i.e.,  datacomputer,  that  is  available  on  the  ARPANET. 
The  datacomputer  will  serve  as  the  actual  on-line  depository  for  much  of  the  raw  and  processed 
seismic  data  generated  by  the  SRO  and  HGLP  programs. 

As  of  this  writing,  not  all  components  of  this  integrated  seismic  system  are  complete.  For 
example,  the  actual  format  descriptions  for  all  the  data  files  which  will  exist  in  the  datacomputer 
are  not  yet  available.  Descriptions  of  the  network  event  processor,  command  and  control  pro- 
cessor, and  the  seismic  information  processor  are  also  not  yet  available.  Apart  from  these 
items,  the  SDIS  contains  a detailed  description  of  all  the  major  components  in  the  SRO  and 
HGLP  programs. 

A directory  has  been  established  on  the  ARPANET  host  computer  located  at  ISIC,  network 
address  number  150,  which  contains  all  the  files  of  the  SDIS.  The  Tenex  protection  for  each 
file  has  been  set  to  allow  any  user  who  is  logged  in  at  ISIC  to  have  read  access  to  these  files. 

This  will  enable  a user  to  read  or  copy  any  file  in  the  directory  (IWWSS)  but  not  allow  him  to 
change  the  files  in  any  way.  Files  within  the  directory  (IWWSS)  can  be  copied  from  the  direc- 
tory to  any  other  host  on  the  ARPANET  by  means  of  the  FTP  programs  located  at  each  site. 
Network  users  cannot  log  into  the  directory  (IWWSS)  but  are  otherwise  free  to  access  the  infor- 
mation stored  there. 
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The  files  within  the  directory  do  vary  in  size  from  a few  pages  to  nearly  4 00  pages.  The 
files  which  reflect  the  current  status  of  the  data  are  short  files  but  will  grow  as  the  seismic 
system  matures.  The  data  status  and  site  status  files  contain  information  that  has  been  taken 
from  a directory  that  is  maintained  by  the  Albuquerque  Seismological  Laboratory.  As  this  di- 
rectory is  updated,  data  will  also  be  incorporated  into  the  (IWWSS)  directory. 

A brief  description  of  each  of  the  major  files  is  given  below.  Persons  with  ARPANET 
access  are  encouraged  to  look  into  the  (IWWSS)  directory.  Persons  without  ARPANET  access 
can  address  a request  to  the  author. 

List  of  Files  Within  the  Directory  (IWWSS) 

Summary- Overview 

List  and  descriptions  of  available  files.  This  list  is  taken  from  the  Summary- 
Overview  file. 

Introduction 

This  file  contains  an  introduction  to  the  documents  available  in  the  directory 
(IWWSS).  Also  covered  is  a brief  description  of  the  components  of  the  Inte- 
grated World  Wide  Seismic  System. 

Index 

A copy  of  the  index  from  each  file  is  contained  here.  It  provides  a good  start- 
ing point  to  determine  what  document  is  appropriate  for  one’s  needs. 

Contacts/People 

People  and  groups  associated  with  the  IWWSS  program.  Indicates  who  to  con- 
tact for  what  kind  of  information. 

Datacomputer -Datafile -Structure 

This  file  gives  proper  data-language  descriptions  of  all  seismic  files  in  the 
datacomputer  and  detailed  English  language  descriptions  of  the  contents  of 
these  files.  When  coupled  with  a proper  users  manual  for  data-language, 
this  file  contains  all  information  needed  by  a user  to  allow  him  to  write  pro- 
grams to  retrieve  data  from  the  datacomputers.  This  file  is  not  yet 
implemented. 

Datacomputer -Version- 1 

This  file  contains  the  programming  manual  for  the  datacomputer.  When  used 
in  combination  with  the  seismic  data  file,  specifications  will  allow  users  to 
obtain  access  to  data  stored  in  the  datacomputer. 

HGLP -Calibrations 

The  most  recent  calibrations  for  the  various  high -gain  long-period  sites  are 
contained  in  this  file.  Older  calibrations  also  will  be  available  in  this  file 
for  historical  purposes. 
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HGLP -General-Description 

General  descriptions  of  the  high-gain  long-period  instrumentation,  seismometer, 
tape  recorders,  and  filters  are  available  in  this  file. 

HGLP-Site-Descriptions 

Specific  site  descriptions  and  installation  reports  are  available  in  this  file.  The 
exact  site  coordinates,  site  geology,  and  descriptions  of  the  installation  condi- 
tions can  be  found  within  this  file. 

SRO  and  HGLP-Data-Status 

For  each  seismic  site,  this  file  lists  when  it  became  operational,  when  major 
changes  in  operation  or  configuration  occurred,  and  what  raw  data  are  available 
from  the  IWWSS  mass  store  and  from  the  original  data  tapes. 

SRO-Calibrations 

The  calibrations  performed  at  the  time  of  installation  and  at  any  time  thereafter 
are  listed  in  this  file. 

SRO-Hardware  and  Software 

The  descriptions  of  the  general  components  of  the  SRO  installations  are  con- 
tained within  this  file.  Those  elements  which  are  common  to  all  SROs  will  be 
described  here.  Items  such  as  the  seismometer,  filters,  computer,  tape 
drives,  etc.,  are  described  in  detail. 

SRO-Site-Descriptions 

The  actual  installation  reports  are  contained  in  this  file.  As  the  report  becomes 
available  it  will  be  included.  Items  covered  include  site  coordinates,  site  geol- 
ogy, and  installation  conditions. 

Seismic -Arrays 

Brief  descriptions  of  the  large  seismic  arrays,  LASA,  NORSAR,  ALP  A,  and 
ILPA  are  contained  in  this  file.  The  actual  seismometer  coordinates  also  can 
be  found  here. 


SIP 

Describes  the  configuration  and  functions  of  the  Seismic  Information  Processor 
(SIP),  its  protocols  and  formats  for  communication  with  the  CCP  and  the  main 
datacomputer;  actual  operational  procedures  including  actions  of  the  data- 
computer,  the  SIP,  the  ARPANET,  the  CCP,  or  the  loss  of  any  data  streams 
into  the  CCP.  This  file  is  not  yet  implemented. 

CCP 

Describes  the  configuration  and  the  Control  and  Communication  Processor;  its 
protocol  and  formats  for  communication  with  the  SIP,  sites  attached  via  the 
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ARPANET,  sites  attached  by  private  lines,  and  the  SDAC  IBM  36  0 computers; 
operational  procedures  including  action  taken  for  all  failure  conditions.  This 
file  is  not  yet  implemented. 


NEP 


Functionally  describes  Network  Event  Processor  and  its  relation  to  the  data- 
computer  and  CCP  in  detail.  This  includes  description  algorithms  and  heu- 
ristics, and  time  scale  and  sequence  for  functions  performed.  Interactions 
between  NEIS,  NEP,  datacomputer,  and  CCP  are  included  in  detail.  This 


file  is  not  yet  implemented. 


R.  M.  Sheppard 


B.  MULTI-AZIMUTH  LARGE  ARRAY  LONG-PERIOD  DATA  BASE 

Several  ongoing  studies  require  a multi-event  data  base  consisting  of  long-period  array 
recordings.  Among  these  seismological  experiments  are:  (1)  a study  to  determine  the  degree 

of  structural  and  mineralogical  velocity  anisotropy  in  the  crust  and  upper  mantle  beneath  ALPA, 
NORSAR,  and  LASA,  (2)  a study  to  determine  the  angular  difference  between  phase  and  group 
wavefronts  for  Love  and  Rayleigh  waves  propagating  across  inhomogeneous  media,  (3)  a study 
to  quantize  the  deviation  in  azimuth  of  surface  waves  from  great  circle  routes  and  to  invert  such 
data  into  path  corrections  for  the  purpose  of  more  uniform  long-period  magnitudes  and  source 
characteristics,  and  (4)  a study  to  compare  the  S-wave  velocity  structure  differences  obtained 
by  the  independent  inversion  of  Love  and  Rayleigh  wave  dispersion  data  at  each  of  the  large 
arrays.  To  these  ends  an  appropriate  data  base  has  been  edited  out  of  the  existing  ASG  seismic 
library,  restored  to  its  original  quality,  and  saved  in  a multifile  FFASTRO  format. 

The  data  format  is  as  follows:  for  each  event,  the  long-period  recordings  of  both  the  ver- 
tical and  horizontal  components  of  each  subarray  triax  are  retained.  The  time  increment  spans 
from  several  minutes  before  the  long-period  P-wave  arrival  to  several  tens  of  minutes  after 
the  high-frequency  tails  of  the  Love  and  Rayleigh  waves.  For  each  of  the  long-period  arrays 
at  LASA,  ALPA,  and  NORSAR,  18  source  events  were  selected  such  that  each  20°  arc  of  azi- 
muth for  each  array  contains  the  great  circle  back  azimuth  of  one  of  the  events.  The  source 
events  for  each  array  were  necessarily  unique  since  the  allowable  depth,  distance,  and  magni- 
tude of  an  event  were  restrained  as  follows:  (1)  the  event  distance  from  the  array  must  be  equal 

to  or  less  than  90°,  (2)  the  event  depth  must  be  less  than  or  equal  to  60  km,  and  (3)  the  event 
magnitude  must  be  equal  to  or  greater  than  5.0  m^.  These  limitations  insure  good  signal-to- 
noise  ratios  over  the  complete  long-period  band.  The  geographic  distribution  for  the  events 
for  the  NORSAR  array  and  the  ALPA  array  are  shown  in  Figs.  VI-1  and  -2,  respectively.  Each 
plot  is  an  equiazimuthal  equidistance  projection  centered  on  the  appropriate  array.  The  LASA 
data  set  is  near  completion.  ^ ^ Needham 

T.  E.  Landers 


C.  SRO  DATA  TAPE  PROCESSING  FACILITY 

A plan  for  PDP-7  Console  expansion  and  other  software  required  to  effectively  handle  SRO 

\ 

data  tapes  was  described  in  a previous  SATS.  This  system  has  been  completely  implemented 
and  is  now  in  operation.  The  last  SATS  discussed  the  first  parts  to  be  completed;  the  final 
programs  that  were  added  to  the  system  were  these: 
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(1)  SROSP,  SROLP,  SR026 

These  programs  run  on  our  PDP-11  and  convert  9-track  SRO  raw  data 
tapes  to  7 -track  tapes  in  the  new  Data  Unit  format  described  in  Ref.  1. 

This  format's  advantage  over  the  SRO  format  is  that  it  includes  more 
descriptive  information  (such  as  instrument  gain  and  orientation)  as 
alphanumeric  data  on  the  tape  itself. 

(2)  Two  new  PDP-7  programs  to  allow  these  Data  Unit  tapes  to  be  used 
as  input  to  the  Analysis  Console  system: 

(a)  DUINIT  (previously  called  >TAB) 

This  program  allows  the  user  to  select  the  data  to  be  viewed  by 
typing  in  an  event  time  and  location.  DUINIT  will  compute  the 
theoretical  arrival  time  at  each  selected  SRO  station  and  will 
set  up  a table  of  arrival  times  to  pass  to  >DUI. 

(b)  >DUI 

>DUI  actually  searches  the  Data  Unit  tape  for  segments  of  data 
that  match  the  stations,  times,  and  orientations  specified  in  the 
table  set  up  by  DUINIT.  When  matching  data  are  found,  >DUI 
reads  them  in  and  initializes  the  Analysis  Console  appropriately. 

As  SRO  tapes  have  been  received  we  have  been  converting  them  to  Data  Unit  tapes  on  the 
PDP-11  and  then  merging  the  tapes  from  various  stations  using  the  PDP-7  DUM  program.  We 
now  have  merged  tapes  available  for  use  starting  with  28  December  1975. 

One  more  software  improvement  is  being  designed:  An  Analysis  Console  plot  program  which 
will  plot  all  the  data  contained  in  an  Analysis  Console  session,  complete  with  alphanumeric  de- 
scriptive information.  (Currently  we  can  only  get  hard  copy  of  the  part  of  the  data  segment  that 
can  be  displayed  on  the  scope  at  any  one  time.) 

L.  J.  Turek 

M.  F.  O'Brien 


D.  PDP-11  SOFTWARE 

An  analysis  of  the  computer  usage  of  Group  22  has  led  to  a reorganization  of  the  PDP-11 
based  computer  systems.  The  analysis  showed  that  the  group  has  a variety  of  data-processing 
requirements.  These  include  interactive  usage  such  as  editing  programs  and  document  prep- 
aration, batch  processing  for  CPU-bound  programs,  both  hard  and  soft  copy  computer  graphics, 
the  maintenance  and  manipulation  of  a large  seismic  data  base,  and  access  to  the  ARPANET. 

All  the  above-mentioned  processing  requirements  are  now  being  performed  on  a number  of 
machines  with  a variety  of  foreign  host  computers  on  the  ARPANET  accessed  through  a USER 
TELNET  running  on  our  PDP-11  systems.  Data  are  transferred  to  and  from  our  site  using  a 
primitive  FTP  program  on  the  PDP-11  system  and  through  the  use  of  the  shuttle  service  be- 
tween Group  22  and  the  Lincoln  IBM  37  0/168. 

The  dispersal  of  these  functions  on  the  various  machines  and  the  quality  of  data  transfer 
connections  between  the  systems  create  a number  of  problems  for  the  user.  Interactive  re- 
sponse time  is  significantly  slower  over  the  ARPANET  than  with  a direct  machine  connection; 
there  is  a lack  of  adequate  network  file  transfer  capability;  the  probability  of  a system  failure 
is  increased  because  the  user  must  go  through  a number  of  computers  to  do  his  processing; 
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the  user  is  faced  with  learning  the  intricacies  of  a number  of  computer  systems;  and,  finally, 
the  local  PDP-11  computers  are  greatly  underutilized.  Because  of  these  problems,  it  was  felt 
that  the  centralization  of  most  of  the  processing  tasks  in  the  Group  22  computers  would  be  de- 
sirable. Only  extremely  large  CPU-bound  programs  and  the  large  seismic  data  bases  would 
remain  on  foreign  host  computers. 

The  development  of  system  software  to  handle  all  these  processing  tasks  would  be  quite 
expensive  and  time  consuming.  It  was  therefore  decided  to  purchase  an  existing  PDP-11  op- 
erating system  and  modify  it  for  our  needs.  The  system  chosen  was  UNIX  developed  by  Bell 

3 

Telephone  Laboratories,  It  was  chosen  because  it  is  an  extremely  well  designed  and  proven 
operating  system  and  contains  .a  large  amount  of  the  required  software. 

With  UNIX  as  a base  it  will  be  possible  to  build  the  other  parts  of  the  desired  system.  It 
will  provide  a software  development  environment  far  superior  to  that  which  we  have  previously 
had.  All  the  standard  features  of  a time-sharing  system,  including  document  preparation  facil- 
ities, already  exist  under  UNIX.  The  ARPANET  software  for  UNIX  is  available  from  the  Uni- 
versity of  Illinois.  An  extremely  sophisticated  interactive  graphics  package  already  exists 
under  UNIX  at  the  University  of  Toronto.  While  this  system  is  more  sophisticated  than  anything 
we  require,  it  can  be  used  as  a basis  for  building  a more  modest  graphics  package. 

A major  problem  which  must  be  solved  is  the  efficient  use  of  computer  resources  to  handle 
the  variety  of  processing  tasks.  It  is  difficult  to  build  a time-sharing  system  which  runs  both 
interactive  and  CPU-bound  programs  efficiently.  Most  systems  which  have  to  run  both  kinds 
of  programs  give  scheduling  priority  to  the  interactive  programs  at  the  expense  of  the  CPU- 
bound  jobs.  A better  solution  for  Group  22  will  be  a dual  CPU  system.  Qie  CPU  will  be  used 
for  all  the  interactive  programs,  while  the  other  will  be  run  in  batch  mode  for  CPU-bound  pro- 
grams. The  machines  will  be  linked  together  through  a data  channel  and  through  a common  disk 
system  (Fig.  VI-3).  The  interactive  processor  will  control  the  system  and  the  user-machine 
interface.  The  batch  CPU  will  run  those  CPU-bound  programs  passed  to  it  by  the  interactive 
processor.  All  data  files  would  be  passed  to  the  batch  program  through  the  common  disk  and 
all  output  will  be  passed  back  to  the  user  through  the  same  disk. 

It  is  felt  that  the  system  described  above  will  significantly  improve  our  computer  user  in- 
terface over  that  which  now  exists.  Users  will  have  a direct  link  to  the  computer,  thereby 
decreasing  response  time  and  increasing  system  reliability.  In  addition,  the  Group  22  re- 
sources will  be  more  efficiently  utilized.  Most  major  hardware  items  are  now  installed.  A 
40-megabyte  disk  system  is  on  order  for  installation  this  summer.  The  UNIX  software  has 
been  ordered  and  will  be  installed  as  soon  as  the  license  to  use  it  is  obtained.  In  the  mean- 
time, a UNIX  system  running  on  another  PDP-11  installation  within  Lincoln  Laboratory  is  being 
used  for  software  development.  T ~ 
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Fig.  VI -1.  Equidistance  projection  for  NORSAR  events. 
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Fig.  VI-2.  Equidistance  projection  for  ALFA  events. 


Fig.  VI-3.  Applied  seismology 
group’s  PDP-11  system  and  its 
relationship  to  the  ARPANET. 
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GLOSSARY 


AIC 

Information  Theoretic  Criterion 

ALPA 

Alaskan  Long  Period  Array- 

ARPA 

Defense  Advanced  Research  Projects  Agency 

ARPANET 

DARPA  Computer  Network 

ASG 

Applied  Seismology  Group 

CAT 

Autoregressive  Transfer  Function  Criterion 

CCP 

Control  and  Communication  Processor 

CPU 

Control  and  Processing  Unit 

EDR 

Earthquake  Data  Reports 

FPE 

Final  Prediction  Error  Criterion 

FTP 

File  Transfer  Protocol 

HGLP 

High  Gain  Long  Period  (Station) 

ILPA 

Iranian  Long  Period  Array 

ILS 

International  Latitude  Service 

ISC 

International  Seismological  Center 

IWWSS 

Integrated  World-Wide  Seismic  System 

LASA 

Large  Aperture  Seismic  Array 

LOD 

Length  of  Day 

LRSM 

Long-Range  Seismic  Measurements 

NEIS 

National  Earthquake  Information  Service 

NEP 

NLS 

Network  Event  Processor 
On-Line  System 

NOAA 

National  Oceanic  and  Atmospheric  Administration 

NORSAR 

Norwegian  Seismic  Array 

NTS 

Nevada  Test  Site 

PDE 

Preliminary  Determination  of  Epicenters 

PEM 

Parametric  Earth  Model 

RMS 

Root  Mean  Square 

SATS 

Semiannual  Technical  Summary 

SDAC 

Seismic  Data  Analysis  Center 

SDIS 

Seismic  Data  and  Information  System 

SIP 

Seismic  Information  Processor 

SPZ 

Short-Period  Vertical  (Seismometer) 

SRO 

Seismic  Research  Observatory 

USCGS 

U.  S.  Coast  and  Geodetic  Survey 

USGS 

U.  S.  Geological  Survey 

WWSSN 

World-Wide  Standardized  Seismograph  Network 
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